Perivascular fluid flow, of cerebrospinal or interstitial fluid in spaces surrounding brain blood vessels, is recognized as a key component underlying brain transport and clearance. An important open question is how and to what extent differences in vessel type or geometry affect perivascular fluid flow and transport. Using computational modelling in both idealized and image-based geometries, we study and compare fluid flow and solute transport in pial (surface) periarterial and perivenous spaces. Our findings demonstrate that differences in geometry between arterial and venous pial perivascular spaces (PVSs) lead to higher net CSF flow, more rapid tracer transport and earlier arrival times of injected tracers in periarterial spaces compared to perivenous spaces. These findings can explain the experimentally observed rapid appearance of tracers around arteries, and the delayed appearance around veins without the need of a circulation through the parenchyma, but rather by direct transport along the PVSs.
Mathematical modelling of ionic electrodiffusion and water movement is emerging as a powerful avenue of investigation to provide new physiological insight into brain homeostasis. However, in order to provide solid answers and resolve controversies, the accuracy of the predictions is essential. Ionic electrodiffusion models typically comprise non-trivial systems of non-linear and highly coupled partial and ordinary differential equations that govern phenomena on disparate time scales. Here, we study numerical challenges related to approximating these systems. We consider a homogenized model for electrodiffusion and osmosis in brain tissue and present and evaluate different associated finite element-based splitting schemes in terms of their numerical properties, including accuracy, convergence, and computational efficiency for both idealized scenarios and for the physiologically relevant setting of cortical spreading depression (CSD). We find that the schemes display optimal convergence rates in space for problems with smooth manufactured solutions. However, the physiological CSD setting is challenging: we find that the accurate computation of CSD wave characteristics (wave speed and wave width) requires a very fine spatial and fine temporal resolution.
What is the right concept of Stokes-Biot stability? The Stokes, Darcy and Biot equations are partial differential equations describing viscous fluid flow, flow in a porous medium, and flow in a poroelastic medium, respectively. Numerical discretizations of the Biot equations are notoriously prone to instabilities and numerical artifacts. In response, Stokes-Biot stability has been introduced as a concept for ensuring Biot stability and convergence even for the challenging cases of low permeabilities and low storage coefficients. The original definition of Stokes-Biot stability relies on both Stokes and Darcy stability. However, ensuring Darcy stability can be highly non-trivial. In this note, we point at how the Darcy stability condition can be relaxed. This allows for a new stability concept: minimal Stokes-Biot stability.
Flow of cerebrospinal fluid (CSF) in perivascular spaces (PVS) is one of the key concepts involved in theories concerning clearance from the brain. Experimental studies have demonstrated both net and oscillatory movement of microspheres in PVS (Mestre et al. (2018), Bedussi et al. (2018)). The oscillatory particle movement has a clear cardiac component, while the mechanisms involved in net movement remain disputed. Using computational fluid dynamics, we computed the CSF velocity and pressure in a PVS surrounding a cerebral artery subject to different forces, representing arterial wall expansion, systemic CSF pressure changes and rigid motions of the artery. The arterial wall expansion generated velocity amplitudes of 60–260 µm/s, which is in the upper range of previously observed values. In the absence of a static pressure gradient, predicted net flow velocities were small (<0.5 µm/s), though reaching up to 7 µm/s for non-physiological PVS lengths. In realistic geometries, a static systemic pressure increase of physiologically plausible magnitude was sufficient to induce net flow velocities of 20–30 µm/s. Moreover, rigid motions of the artery added to the complexity of flow patterns in the PVS. Our study demonstrates that the combination of arterial wall expansion, rigid motions and a static CSF pressure gradient generates net and oscillatory PVS flow, quantitatively comparable with experimental findings. The static CSF pressure gradient required for net flow is small, suggesting that its origin is yet to be determined.
In these studies, we consider systems of equations describing flow in multiple, interacting fluid networks with different permeabilities. The resulting discrete systems of equations can easily become very large, and efficient solution techniques are key. In cases where the inter-network interactions dominate the intra-network permeabilities, we observed that the systems of equations became increasingly hard to solve. In particular, the number of solver iterations would increase with increasing discrete resolution and interaction parameters. However, by creating an algorithm for defining tailored changes of variables, we were able to define robust preconditioners with uniformly bounded iteration counts.
Personal favorite: The lovely trick for how to define parameter-tailored changes of variables via classical matrix analysis including simultaneous diagonalization by congruence.
Mathematical models in biology involve many parameters that are uncertain or in some cases unknown. Over the last years, increased computing power has expanded the complexity and increased the number of degrees of freedom of many such models. For this reason, efficient uncertainty quantification algorithms are now needed to explore the often large parameter space of a given model. Advanced Monte Carlo methods such as quasi Monte Carlo (QMC) and multilevel Monte Carlo (MLMC) have become very popular in the mathematical, engineering, and financial literature for the quantification of uncertainty in model predictions. However, applying these methods to physiologically relevant simulations is a difficult task given the typical complexity of the models and geometries involved. In this paper, we design and apply QMC and MLMC methods to quantify uncertainty in a convection-diffusion model for tracer transport within the brain. We show that QMC outperforms standard Monte Carlo simulations when the number of random inputs is small. MLMC considerably outperforms both QMC and standard Monte Carlo methods and should therefore be preferred for brain transport models.
Mathematical models for excitable cells are commonly based on cable theory, which considers a homogenized domain and spatially constant ionic concentrations. Although such models provide valuable insight, the effect of altered ion concentrations or detailed cell morphology on the electrical potentials cannot be captured. In this paper, we discuss an alternative approach to detailed modelling of electrodiffusion in neural tissue. The mathematical model describes the distribution and evolution of ion concentrations in a geometrically-explicit representation of the intra- and extracellular domains. As a combination of the electroneutral Kirchhoff-Nernst-Planck (KNP) model and the Extracellular-Membrane-Intracellular (EMI) framework, we refer to this model as the KNP-EMI model. Here, we introduce and numerically evaluate a new, finite element-based numerical scheme for the KNP-EMI model, capable of efficiently and flexibly handling geometries of arbitrary dimension and arbitrary polynomial degree. Moreover, we compare the electrical potentials predicted by the KNP-EMI and EMI models. Finally, we study ephaptic coupling induced in an unmyelinated axon bundle and demonstrate how the KNP-EMI framework can give new insights in this setting.
Mixed dimensional partial differential equations (PDEs) are equations coupling unknown fields defined over domains of differing topological dimension. Such equations naturally arise in a wide range of scientific fields including geology, physiology, biology and fracture mechanics. Mixed dimensional PDEs are also commonly encountered when imposing non-standard conditions over a subspace of lower dimension e.g. through a Lagrange multiplier. In this paper, we present general abstractions and algorithms for finite element discretizations of mixed domain and mixed dimensional PDEs of co-dimension up to one (i.e. nD-mD with |n-m| <= 1). We introduce high level mathematical software abstractions together with lower level algorithms for expressing and efficiently solving such coupled systems. The concepts introduced here have also been implemented in the context of the FEniCS finite element software. We illustrate the new features through a range of examples, including a constrained Poisson problem, a set of Stokes-type flow models and a model for ionic electrodiffusion.
Uncertainty in diffusion parameters substantially impact the amount of tracer in gray and white matter, and the average tracer concentration in gray and white subregions.
Even with an uncertainty in the diffusion coefficient of a factor three, and a resulting fourfold variation in white matter tracer enhancement, discrepancies between simulations of diffusion and experimental data are too large to be attributed to uncertainties in the diffusion coefficient alone.
A convective velocity field modelling the glymphatic theory did not increase tracer enhancement in the brain parenchyma compared to pure diffusion. However, when a large-scale directional structure was added to this glymphatic velocity field, tracer inflow increased.
Infusion testing is a common procedure to determine whether shunting will be beneficial in patients with normal pressure hydrocephalus. The method has a well-developed theoretical foundation and corresponding mathematical models that describe the CSF circulation from the choroid plexus to the arachnoid granulations. Here, we investigate to what extent the proposed glymphatic or paravascular pathway (or similar pathways) modifies the results of the traditional mathematical models.
We used a two-compartment model consisting of the subarachnoid space and the paravascular spaces. For the arachnoid granulations, the cribriform plate, capillaries and paravascular spaces, resistances were calculated and used to estimate flow before and during an infusion test. Next, pressure in the subarachnoid space and paravascular spaces were computed. Finally, different variations to the model were tested to evaluate the sensitivity of selected parameters.
At baseline, we found a very small paravascular flow directed into the subarachnoid space, while 60% of the fluid left through the arachnoid granulations and 40% left through the cribriform plate. However, during the infusion, paravascular flow reversed and 25% of the fluid left through these spaces, while 60% went through the arachnoid granulations and only 15% through the cribriform plate.
The relative distribution of CSF flow to different clearance pathways depends on intracranial pressure (ICP), with the arachnoid granulations as the main contributor to outflow. As such, ICP increase is an important factor that should be addressed when determining the pathways of injected substances in the subarachnoid space.