Article Information

• PDF (642 kb)
• Full Text with Large Figures
• Cited by in Scopus (2)

PubMed

• Articles by Ravi K. Nandigam
• Articles by Daniel M. Kroll

Related Articles

• The Invagination of Excess Surface Area by Shrinking Neurons...

  The Invagination of Excess Surface Area by Shrinking Neurons Biophysical Journal, Volume 85, Issue 1, 1 July 2003, Pages 223-235 C.E. Morris, J.A. Wang and V.S. Markin Abstract Over most of their surface, neurons are surrounded by a narrow extracellular gap across which they make adhesive cell-cell contacts. Thus constrained, how do they regulate their geometry when osmotically perturbed? Specifically, are there any interesting consequences of local osmosis in such conditions? Using confocal imaging of shrinking neurons in culture, we observe water exiting into the cell-substratum gap. This water efflux generates a hydrostatic pressure that, at discrete (low adhesion) sites, causes the neuron's excess plasma membrane to invaginate, thus compensating for shrinkage with a pseudo-intracellular volume. To identify the minimal requirements of the process, a compartment/flux model was constructed. It comprises, essentially, a large liposome adhering in a labyrinthine fashion to a substratum. The model predicts that invaginations form at the cell-substratum interface under the influence of local osmosis, provided that adhesion across the gap is neither too tight nor too loose. Local osmosis in the central nervous system, in contrast to epithelia, is usually considered a mishap, not a physiological opportunity. We postulate, however, that local osmotic forces acting in conjunction with confined extracellular spaces could be harnessed in service of surface area, shape, and volume regulation when intense neural activity alters a neuron's osmotic balance. Abstract | Full Text | PDF (679 kb)

• Extracellular space structure revealed by diffusion analysis...

  Extracellular space structure revealed by diffusion analysis Trends in Neurosciences, Volume 21, Issue 5, 1 May 1998, Pages 207-215 Charles Nicholson and Eva Syková Abstract The structure of brain extracellular space resembles foam. Diffusing molecules execute random movements that cause their collision with membranes and affect their concentration distribution. By measuring this distribution, the volume fraction () and the tortuosity () can be estimated. The volume fraction indicates the relative amount of extracellular space and tortuosity is a measure of hindrance of cellular obstructions. Diffusion measurements with molecules <500 show that ≈0.2 and ≈1.6, although some brain regions are anisotropic. Molecules ≥3000 show more hindrance, but molecules of 70000 can move through the extracellular space. During stimulation, and in pathophysiological states, and change, for example in severe ischemia = 0.04 and = 2.2. These data support the feasibility of extrasynaptic or volume transmission in the extracellular space. Abstract | Full Text | PDF (2576 kb)

• Poly[N-(2-hydroxypropyl)methacrylamide] Polymers Diffuse in ...

  Poly[N-(2-hydroxypropyl)methacrylamide] Polymers Diffuse in Brain Extracellular Space with Same Tortuosity as Small Molecules Biophysical Journal, Volume 80, Issue 1, 1 January 2001, Pages 542-548 Š. Prokopová-Kubinová, L. Vargová, L. Tao, K. Ulbrich, V. Šubr, E. Syková and C. Nicholson Abstract Integrative optical imaging was used to show that long-chain synthetic poly[-(2-hydroxypropyl)methacrylamide] (PHPMA) polymers in a range of molecular weights from 7.8 to 1057kDa were able to diffuse through the extracellular space in rat neocortical slices. Tortuosity (square root of ratio of diffusion coefficient in aqueous medium to that in brain) measured with such polymers averaged 1.57, a value similar to that obtained previously with tetramethylammonium, a small cation. When PHPMA was conjugated with bovine serum albumin (BSA) to make a bulky polymer with molecular weight 176kDa, the tortuosity rose to 2.27, a value similar to that obtained previously with BSA alone and with 70-kDa dextran. The method of image analysis was justified with diffusion models involving spherical and nonspherical initial distributions of the molecules. Abstract | Full Text | PDF (387 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 10, 3368-3378, 15 May 2007

doi:10.1529/biophysj.106.095547

Biophysical Theory and Modeling

Three-Dimensional Modeling of the Brain's ECS by Minimum Configurational Energy Packing of Fluid Vesicles

Ravi K. Nandigam^* and Daniel M. Kroll^†, ,  

^* School of Chemical Engineering, Purdue University, West Lafayette, Indiana
^† Department of Physics, North Dakota State University, Fargo, North Dakota

Address reprint requests to D. M. Kroll, Tel.: 701-231-8968 or 231-8974.

The interstitial space between cells in the brain is called the extracellular space (ECS) 1. Diffusion in the ECS is an essential component in many processes ranging from the delivery of glucose to cells to novel intracellular communication; it is also the primary mechanism for drug delivery to the brain and is important for measurements such as diffusion-weighted magnetic resonance imaging.

The brain is composed of a heterogeneous packing of a large number of different cell types. While the cell bodies of brain cells typically have a diameter in the range of 5–50μm, many of the cellular elements in the central nervous system are smaller. A typical brain cell is enclosed in a lipid bilayer membrane containing various receptors and protein channels, and there are numerous extensions from the cell body—called processes—which play a fundamental role in intercellular communications 2. Although these cells are densely packed, there is a small ECS between each cell. The cells do not appear to be in direct contact, except at gap junctions, and it is generally estimated that the width of the gap between the cells is on the order of 200Å 3. There is a large amount of experimental data that clearly shows that the ECS comprises ∼20% of the total volume of the brain 1. The ECS is filled with an electrolyte containing various ions and a small number of larger organic molecules. The collection of long-chain molecules in the ECS forms what is called extracellular matrix 4. Many of these molecules are tethered to the cell membranes, but others float freely in the ECS. Although there is considerable interest in the extracellular matrix, its density and composition, as well as its influence on the diffusive processes in the ECS, are not well known.

Most of what is known about the structure of the ECS comes from measuring the diffusion constants of various molecules. This data can be summarized in terms of the behavior of the tortuosity, λ, defined here as the square-root of the ratio of the diffusion constant of the molecule in a free medium (water or dilute gel), D, over the value in the ECS (D*), i.e., In healthy brain tissue, λ is typically ∼1.6, but can be as large as 1.9–2 when there are pathologies that involve cellular swelling.

The extra cellular space has a very complex structure, since the shapes of the constitutive cells are complex and the connecting passages are tortuous and have a random connectivity. There are two distinct factors contributing to the measured tortuosity. The first is purely geometrical: in the ECS, the diffusing particle has only a limited number of tortuous paths to follow, requiring a longer time to diffuse a given distance 5. This contribution to the tortuosity is called the geometric tortuosity. Another possible contribution comes from an enhanced viscosity in the ECS (called constitutive effects) caused by the impeded movement of the diffusing molecule by the extracellular matrix 6,7. Only the geometric dependence of tortuosity will be considered in this article.

There are also a considerable amount of experimental data on the behavior of the volume fraction of the ECS, α, and the tortuosity, λ, in osmotically stressed brain tissue. In recent experiments, the size of the ECS was controlled by varying the osmotic pressure of media bathing isolated brain tissue 8,9. Measurements on rat cortex have shown 10 that while α decreases under hypo-osmolarity and increases with hyperosmolarity, λ increases in hypotonic media, but reaches a plateau after first decreasing slightly in hypertonic media, in contradiction to the usual expectation that λ should monotonously decrease with increasing α.

There have been a large number of attempts to compute the tortuosity from models representing the brain structure. Rusakov and Kullmann 6 constructed a computational model based on differential geometry to determine the tortuosity of brain tissue. Their model involved representing the ECS as a three-dimensional random porous structure and computing the increase in the mean local path length of the diffusing molecule in the tissue as compared to the free medium. However, as Nicholson 1 pointed out, such geometric approaches suffer from the drawback that it is hard to correctly weight the contributions of different paths.

Chen and Nicholson 11 used the area outside regular arrays of convex impermeable surfaces with rounded corners to model the ECS in two dimensions. The tortuosity for these models was computed by transforming the problem into one that involved solving a simple Laplace equation in the interstitial region, by applying principles of homogenization theory 12. Using these models, they studied the dependence of the tortuosity on the separation 2h between the planar surfaces in the arrays and the radius of curvature a of the corners (see Fig. 1). They also considered random arrangements of cells of random shape in two dimensions. On the basis of their results, they argued that the behavior of λ as the brain underwent osmotic stress could be due to fact that the brain cells change their shape with changes in osmolarity. They argued that it is possible that λ can remain constant with increasing α if there are “lakes” at the junctions between the cells, which locally trap diffusing molecules, thereby offsetting the increased diffusion in the expanding free space.

Display large version of this figure

Figure 1

(A) Two-dimensional regular arrangement of square cells with rounded corners. (B) Detail of the periodic repeat unit used for testing the random walk algorithm and various implementations of zero flux boundary conditions. The figure has been redrawn from Chen and Nicholson 11.

More recently, Tao and Nicholson 13 considered several model geometries of periodic arrays of convex polytopes in three dimensions. In these models, the ECS was an interconnected planar array of uniform thickness. For these geometries, they obtained the simple relationship which is the same as a result Maxwell derived for a dilute suspensions of spheres 14, and argued that three-dimensional models constructed from convex polytopes do not capture the critical structural features of the brain's ECS. In subsequent work, which incorporated concave elements into the model 15,16,17,18, it was found that geometrical tortuosities very similar to those measured in experiments were obtained if the “dead space” created by the concave invaginations were sufficiently large. Hrabětová et al. 19 presents some experimental justification for dead-space microdomains, as well as some suggestions about how they may form.

In this article, we model the ECS network in three dimensions by relaxing an array of fluid membrane vesicles. While this first-principles approach yields only a crude approximation of the ECS of brain tissue, two-dimensional cross-sections of the resulting networks look very similar to the random networks studied in Chen and Nicholson 11, and it is likely that shape changes at ECS network junctions induced by osmotic challenge can be modeled reasonably accurately. Since the initial arrays of spherical vesicles have porosities in the range 32–40%, depending on the degree of polydispersity 20, we expect—and do find—that the void space at junctions in the resulting ECS networks is relatively large. The ECS networks we generate consist of curved passages along multicell contact lines connecting these junctions. Diffusivities in this network are measured using the particle tracking algorithm developed to study dispersion in beadpacks 21,22. The influence of osmotic stress is studied by imposing a pressure difference between the interior and exterior regions of the vesicles.

The brain is composed of a complex array of cells with a wide range of size, structure, and function. As a result, the connecting passages in the ECS are tortuous, the connectivity is random, and there are many relevant length scales. Furthermore, many of the cellular elements are not spherical, but have a tubular or sheetlike structure. A first-principles modeling of diffusion in the brain's ECS is therefore extremely difficult, and the complexity of the problem is the primary reason that most efforts to model the ECS have utilized regular arrays of convex polytopes.

Cells are bounded by a lipid bilayer membrane, and the interior of the cell contains the cytoskeleton, a complex filamentous network comprised of microtubules, actin filaments, and intermediate filaments. Together, these structures provide shape and mechanical integrity for the cell. While the elastic properties of individual cells are not well understood, the mechanical properties of lipid-bilayer membranes have been studied intensively recently, and it has been shown that their elastic and thermal properties can be well described using a theory based on the curvature elasticity of thin shells. For example, the characteristic discocyte shape of a red blood cell can be obtained by minimizing the curvature energy with the constraints of fixed membrane area and a particular fixed enclosed volume.

Ideally, if we knew how to model the mechanical properties of these cells, we could relax random arrays of the cells to form structurally representative examples of the brain and its extracellular space. Various techniques such as homogenization theory or random walk algorithms could then be used to determine the tortuosity as a function of the porosity. However, as noted above, the structure and elastic behavior of cells is quite complicated, and it is not clear how it should be modeled. What we have done, therefore, is to use a simple model of the curvature elasticity of thin shells to describe the mechanical properties of the cells. While this is, at best, only a very crude approximation to the elastic properties of real cells, it enables us to study the dependence of the tortuosity on pore volume fraction and osmotic stress in both regular and random arrays of vesicles with well-defined elastic properties. In this article, we do not consider the influence of the extracellular matrix. Although measurements of the molecular weight dependence of the tortuosity indicate that the extracellular matrix may significantly affect the values of the tortuosity in certain cases, due to the limited understanding of viscous effects on tortuosity, we do not consider them here. We hope that the current approach will provide insight into the influence of various structural factors on the geometric tortuosity, and serve as a guide for more detailed future studies.

In the current approximation, the shape of the cell is determined by the curvature elasticity of the lipid bilayer surrounding the cell; the molecular properties, architecture, and the interactions of the membrane constituents enter only through the functional form of the elastic energy and the values of the elastic moduli.

A lipid membrane is free to adjust its surface area so as to minimize the free energy of the amphiphiles. On the timescale of a typical experiment, the number of lipids in a membrane does not change, so that the area of the membrane is constant and the interfacial tension is zero in the unstressed state. The shape and fluctuation spectra of lipid membranes are determined by its curvature elasticity.

As already noted, the Hamiltonian of fluid membranes is not only invariant under rotations and translations, but also under reparameterizations. This additional invariance is due to the fluid structure, which does not allow a preferred coordinate system, and therefore cannot support shear stress. Fluid membranes are compressible, but the compressibility modulus is usually rather large, so that they are often studied in the incompressible limit, where the membrane area is fixed. In this case, the only contribution to the configurational energy is the bending energy 23,24,25

(1)

We model a membrane as a triangulated network of particles that form a regular two-dimensional array embedded in d=3 spatial dimensions. Vesicles are closed sheets with the topology of a sphere.

We employ a simple string-and-bead model 27 in which the particles or vertices of the triangular network interact via the tethering potential

(2)

The potential V(r) acts only between tethered nearest neighbors; it ensures that the distance between nearest neighbors is In simulations, self-avoidance can be guaranteed by placing a particle at each vertex that is sufficiently large that it cannot pass through the network mesh. Because of the large values for the bending rigidity we employ, self-intersection of individual vesicles does not occur, and this was not necessary. Nevertheless, it is convenient to include the hard sphere potential

(3)

To model fluid membranes, this network model has to be modified to allow for the diffusion of vertices in the membrane. This is done by making the connectivity of the network a dynamic variable. The simplest way to do so is to cut and reattach tethers between the four beads of two neighboring triangles 28,29,30,31,32,33. To maintain the triangular nature of the network, a bond-flip is only allowed if the initially connected vertices have at least four neighbors each. Also, a bond-flip is only possible if the distance of the two initially disconnected vertices falls within the tether length. The bond-flip algorithm has the advantage that it preserves both the two-dimensional connectivity and the topology of the network 27.

We assume that the spontaneous curvature of the membrane is zero, and use a discretization of the squared Laplacian form of the bending energy

(4)

A general introduction to methods for discretizing operators on triangulated random surfaces is given by Itzykson 34. On a triangulated surface, the mean curvature at node i is

(5)

(6)

Since n‖ΔR for surfaces embedded in three dimensions, Eq. (5) implies that the Laplacian squared bending energy can be written as 34,35

(7)

Using the discrete free energy Eq. (7), we have generated various three-dimensional arrays of nonintersecting vesicles by minimizing the configurational free energy. A general introduction to the Monte Carlo procedure we employed can be found in Gompper and Kroll 27. Both regular and random networks have been generated. In the first case, a Monte Carlo simulation of a single vesicle was performed in a cubic box. The size of the box was periodically allowed to shrink, while making sure that the vesicle does not intersect the box surface. This was done by checking that none of the tethers intersect the box surfaces. The value κ/k_BT=30 was used for the bending rigidity, ensuring that that vesicle surface is smooth. This procedure was continued until a porosity of 20% was attained. Simulations performed using larger values of the bending rigidity yielded configurations consistent with those obtained for this value of κ. A simple cubic periodic array was then obtained by replicating this elementary cell. Osmotically stressed periodic arrays of various porosities were generated in the same way while applying a pressure difference between the vesicle interior and exterior, keeping the box size fixed.

Random networks were generated using a similar procedure after starting with a random close-packed array of spherical vesicles in a simulation cell with periodic boundary conditions. As the size of the simulation cell was decreased, the vesicles deformed, creating a tightly packed random array of highly deformed vesicles. During this procedure, self-avoidance was ensured by making sure that a vesicle's tethers do not intersect the surface triangles of adjacent vesicles. Arrays with a wide range of porosities were generated in this way. We did not consider osmotically stressed random arrays.

The interstitial spaces separating the vesicles were then used to model the ECS.

Lagrangian particle tracking methods are more efficient than finite element methods for solving the convection-diffusion equation in problems where Peclet numbers are very large and/or the geometries are complex 36. Particle tracking methods are more stable, easy to implement, and free from numerical dispersion and mesh generation problems. In the current application, the Peclet number is zero. However, the three-dimensional geometries we consider are very complex, and mesh generation and refinement problems prevented us from using finite element methods in three dimensions. We therefore used a random walk method, which is the most straightforward implementation of the Lagrangian particle tracking method. In the random walk model we use, the position of a tracer particle is determined by the following stochastic differential equation 21,

(8)

(9)

It is well known that random walk algorithms describe a diffusion process. However, there is no consensus regarding the correct implementation of certain boundary conditions for particle-based methods. In this work, we use both periodic and zero-flux boundary conditions. Although the implementation of periodic boundary conditions is straightforward, several procedures have been used in the literature to describe zero-flux boundary conditions 21. The appropriate choice is analyzed in the following.

Proposed implementations of the zero-flux boundary condition are

1. Specular reflection. The trajectory, r(t) of a tracer particle undergoing specular reflection at the x=0 plane is given by

(10)

where R is the mirror reflection operator and Ω is the diffusion domain. The value is the unit vector normal to the x=0 plane 36.
2. Rejection. The move is rejected if r(t)+ΔW(t)∉Ω, i.e., for that time step Δt the particle remains at its previous position 37,38.
3. Multiple rejection. The increment ΔW(t) is rejected and a new increment is found until r(t)+ΔW(t) ∈ Ω 39.
4. Interruption. The particle stops at the boundary and its clock is incremented by γΔt with γ given by

(11)

It is shown in Szymczak and Ladd 36 that some of these schemes (multiple rejection and interruption) fail to impose the zero-flux boundary condition correctly. For the current application, we determine the most accurate procedure by performing random walk simulations in a well-defined two-dimensional geometry whose tortuosity can be accurately determined by an alternate method. The procedure that yields the most accurate results will be used in subsequent three-dimensional simulations.

We have chosen to use the two-dimensional geometry shown in Fig. 1, which is one of the lattice arrangements of cells studied by Chen et al. 11. Accurate estimates for the tortuosity have been determined by applying the homogenization method discussed in detail in Chen and Nicholson 11. Since the geometry is simpler in this case, mesh generation problems did not occur. The construction of the two-dimensional geometry and the mesh generation is done in GAMBIT (Fluent Inc., Lebanon, NH), and the tortuosity is computed by numerically solving the final Laplace equation using FLUENT (Fluent, Inc.) software. FLUENT's user-defined scalar feature has been used for defining the variable ω of Chen and Nicholson 11, and the default parameters in FLUENT—velocity, viscosity, etc.—are set to zero. Also, the default criterion for convergence in FLUENT has been used.

Results for the tortuosity as a function of a, the radius of curvature of the corners, for fixed a/h=1 and 1.5, are shown in Fig. 2. For fixed a/h, it can be seen that λ(a) is linear in a. Although the slope of the straight line depends on the value of the ratio a/h, the intercept with λ-axis is the same, e.g., 1.414, in both cases. For fixed a/h, the limit a → 0 results in a network of tubes, which is known to have the tortuosity 43.

Display large version of this figure

Figure 2

Plots of λ vs. a for a/h=1.5 and a/h=1. The value a is the curvature radius of the cell corners in the two-dimensional array shown in Fig. 1. Tortuosities were computed using the approach based on homogenization theory.

Four sets of random walk simulations were performed using the same geometry, with each set corresponding to one of the four implementations of the zero flux boundary conditions discussed above. The region lying within the unit square and outside the shaded region is the interstitial domain occupied by the diffusing particles. Periodic boundary conditions are imposed at the edges of the unit cell (to model an infinite periodic system), and no-flux boundary conditions were implemented on the solid boundaries. At the start of the simulation, the domain is uniformly filled with large number of tracer particles. At each time step, the coordinates of the tracer particles are updated according to

(12)

Display large version of this figure

Figure 3

Repeat unit used in the random walk simulations of diffusion in the regular array shown in Fig. 1. The figure also illustrates the implementation of periodic boundary conditions for two hypothetical tracer particles.

The second moment of the relative global displacement of the particles, which is needed to determine the effective diffusion constant, is given by

(13)

(14)

Results for the tortuosities for a/h=1 are plotted in Fig. 4, where it can be seen that only the specular reflection implementation of the zero-flux boundary condition yields the correct results for the tortuosities for all values of a. The expected linear relation between λ and a is only observed for this case. Table 1 contains a comparison of results for λ(a) obtained using the homogenization method and the random walk algorithm using specular reflection boundary conditions.

Display large version of this figure

Figure 4

Plots of λ vs. a. Values of the tortuosity were obtained by performing random walk simulations of diffusion in the two-dimensional lattice arrangement of Fig. 1 for various implementations of the zero-flux boundary condition at the cell surface. The value a/h=1, where a is the radius of curvature of the cell corners and h is half the gap width between the cells.

As can be seen in Fig. 4, for other implementations of the no-flux boundary condition, results for the tortuosity begin to deviate from the correct result as a approaches zero, i.e., as the tubes become narrower. For the rejection boundary condition, the move is rejected when a particle hits the solid surface. Particles near solid surfaces therefore move slower than those in the bulk, and the tortuosity is somewhat larger for all tube diameters. As the tubes become narrower and narrower, more and more of the particles hit a solid boundary during an update, and an ever increasing number of moves are rejected. Consequently, particles in these regions tend to remain there, and the tortuosity increases. In the case of the multiple rejection boundary condition, new increments are attempted until one is found that does not intersect the surface. In the narrow tubular regions, this means that the diffusion becomes one-dimensional as h goes to zero. In addition, since the particles cannot be deflected by the walls, they tend to continue to move along their line of motion at the junctions in this limit. The diffusion process therefore becomes one-dimensional as h → 0, and the tortuosity approaches 1 ( since d=1 here). Finally, for the interruption boundary condition, the overall move is a combination of two or more random variables; the variance of the particle's displacement (and hence overall displacement of the particles) is therefore increased, and the computed tortuosity is smaller than the correct value.

The convergence properties of the Euler approximation for linear stochastic equations depend on whether one seeks to approximate the moments of the particle distribution, or to satisfy an error criterion for particle position 44. The former requirement is known as weak convergence, and the rate is linear, both in the number of particles, N, and in the time step, Δt; the latter is known as strong convergence, and the rate is 21. In the present work, since we are interested in finding tortuosities in different porous media, which would require computing the second moment of the particle displacements, the weak convergence criterion is used.

The simulations in two dimensions were performed using between 10⁴ and 2×10⁴ tracer particles. Generally, it was found that 10,000 particles was sufficient. The size of time step (Δt) was chosen so that a further reduction in Δt resulted in no changes in the measured tortuosity. The size of the time step size depended on the value of a. The appropriate choice of Δt decreased from 4×10⁻⁴ to 5×10⁻⁵ as a was reduced from 0.3 to 0.005. Because of the initial uniform distribution of tracer particles, converged results were obtained using no more than 1000 time steps. This was sufficient to produce results for which scaled linearly in time.

Three-dimensional simulations were performed in the void spaces of both the regular and random vesicle packs generated using the procedure described above. Figure 5 and Figure 6 show the elementary repeat unit for the cubic array and the random network. In both cases, the unit cell is a cube with edges of length 1. At each time step, the coordinates of the tracer particles are updated according to

(15)

Display large version of this figure

Figure 5

Single relaxed vesicle in a periodic unit cube. The unit cube is replicated periodically to form the network. Neighboring cells are in contact along their flat surfaces. Diffusion occurs in the void space external to the cubes.

Display large version of this figure

Figure 6

A random array of 18 highly deformed vesicles forming a random porous network. The unit cell is the unit cube, which is replicated periodically to form the network. Diffusion occurs in the void space between the vesicles. (The number of vesicles in the figure is >18 because periodic images of some of the vesicles were included to avoid cuts through the vesicles.)

At the start of the simulations, the pore space was populated with a random uniform distribution of tracer particles. The number of tracer particles was typically on the order of 50,000. At each time step, the coordinates of the tracer particles were updated according to Eq. (15) and specular reflection boundary conditions were used if the particle intersected a solid surface.

To improve the efficiency of the algorithm when determining whether a particle's trajectory intersects a vesicle surface during a time step, the unit simulation cell is divided into several subcells. Each subcell has a linked list of all surface triangles located in the subcell 46. When a particle's coordinates are updated, it is then only necessary to check for intersections with surface triangles in the same and nearest-neighbor subcells. This reduces the computational time by a factor proportional to the number of subunit cubes.

The number of time steps in the simulation is chosen sufficiently large that the slope of a plot of vs. t reaches its asymptotic value. For the configurations considered in this article, it was found that 1000 time steps were sufficient for simulations involving 50,000 tracer particles. The time increment, Δt, was initially taken to be 4×10⁻⁵, so that the particle's displacement during a time step was ∼5% of the size of the unit simulation cell. In a given simulation, Δt was gradually reduced until the asymptotic slope of a plot of vs. t did not change if the time step is further reduced.

The pore volume fraction, α, does not characterize the structure of the porous domain. Additional information is contained in the pore-size distribution, P(s), which is the probability that a small volume element of pore space is located a (nearest) distance s from a pore surface 47. To determine P(s), a large number (∼100,000) of uniformly distributed points are randomly chosen in the void space and the closest distance s of each point to any pore surface is measured. P(s) is the normalized sum of the number of points, which lie at a distance s±ds from the solid surface. Since the normalized distribution satisfies P(s) has the dimension of inverse length. P(0) is the ratio of the pore surface area to the pore volume, S/V_p≡1/R_h, where R_h is called the mean hydraulic radius, a quantity that can be measured in diffusion nuclear magnetic resonance experiments 48.

Results for the pore-size distribution in both the periodic and random porous networks we generated are shown in Figure 7 and Figure 8, respectively. When comparing the results, it is important to remember that in both cases, the periodic simulation cell has a linear dimension of 1, and s is measured in these units. The sharp peak of P(s) at s very close to 0, indicates that the cells are in very close contact with each other. This remains true even as the pore volume is increased; however, the area of contact is significantly reduced in this case, implying a nonuniform shrinkage of cells—as proposed in Chen and Nicholson 11 during osmotic stress.

Display large version of this figure

Figure 7

Pore-size distributions of the regular porous network for a range of pore volume fractions. P(s) is the fraction of the porous region that is distance s from the nearest cell wall.

Display large version of this figure

Figure 8

Pore-size distributions of the random porous network for a range of pore volume fractions. P(s) is the fraction of the porous region that is distance s from the nearest cell wall.

Our results for the geometric tortuosities as a function of the pore volume fraction for both the regular and random packings are shown in Fig. 9. The range of pore volume fractions we considered is 0.11–0.27 for regular networks and 0.10–0.30 for the random packing. First, note that the tortuosities we calculate for the cubic array are significantly larger than those obtained by Tao and Nicholson 13 in their study of uniformly spaced convex cells. The reason for this is that the networks considered in Tao and Nicholson 13 consisted of polytopes with sharp corners. In contrast, the configurations we obtain by minimizing the curvature energy (see Fig. 5) have rounded vertices. In fact, inspection of Fig. 5 shows that the curvature at the cube vertices is smaller than along the edges; as a consequence, the pore space at the junctions is comparatively large.

Display large version of this figure

Figure 9

Tortuosities obtained from random walk simulations in the regular and random porous networks as a function of the pore volume fractions.

The effect of having junctions with a larger radius of curvature than the typical dimension of the connecting tubular regions can be seen in Fig. 10, which contains a plot of the tortuosity of the two-dimensional network shown in Fig. 1 as a function of h, for fixed a=0.4. For h>0.1, the tortuosities plotted in Fig. 10 are slightly smaller than those shown in Fig. 2 for h=a. For smaller h, however, the tortuosities increase dramatically, since diffusing particles spend an ever-increasing amount of time trapped at the junctions as h → 0. The ratio of the characteristic size of junctions to the connecting tubular regions is therefore a primary factor that determines the tortuosity, and explains the differences between our results and those of Tao and Nicholson 13.

Display large version of this figure

Figure 10

Plot of λ as a function of h for fixed a=0.4. a is the radius of curvature of the cell corners in the regular two-dimensional array shown in Figure 1h is half the gap width between the cells. Tortuosities were computed using the approach based on homogenization theory.

The results presented in Fig. 9 also show that the random network has tortuosities that are smaller than those for the regular array. This is probably due to the fact (see Fig. 11) that while the cubic array has an organized structure consisting of narrow channels and junctions, the random array does not (at least at the current resolution). Indeed, several regions of the void space in the random pack are surprisingly large; higher resolution studies (in which the number of surface triangles in a vesicle is increased) will be required to determine whether this remains true in the continuum limit. Insufficient resolution could also be a reason that the measured tortuosities for the random pack are lower than observed in experiment. The size of a typical surface triangle is 0.05 for cubic regular array (Fig. 5), and 0.1 for random pack (Fig. 6). This means that length scales smaller than this size are not resolved correctly. While this does not seem to be a problem for the cubic array, it is probable that higher resolution studies would show that the channels connecting the larger junctions in the random pack are much smaller than predicted here. On the other hand, the larger junctions are reasonably well resolved. Since an overestimation of the size of small channels reduces the tortuosity, we expect the tortuosity to increase for higher resolutions.

Display large version of this figure

Figure 11

Comparison of sample cross-sections of the random porous network (left) and regular network (right). In both cases, the pore volume fraction is ∼0.2.

It can also be seen that for the random network, the magnitude of the slope of the tortuosity versus pore volume fraction plot decreases with increasing α. This is what is observed in experiment. In fact, an overlay of the tortuosities of the random network with the experimental data (plotted in Fig. 12) shows that the random network does indeed provide a reasonably accurate representation of the pore volume fraction dependence of the tortuosity. It is important to note, however, that for pore volume fractions larger than the values studied, the vesicles are essentially spherical, and there is no plateau. For the plateau to extend to much larger values of porosity—as observed in experiment—it would probably be necessary to consider packings of polydisperse vesicles with nonspherical rest shapes.

Display large version of this figure

Figure 12

A comparison of simulation and experimental results for the tortuosity as a function of the pore volume fraction, α. The simulation results (right vertical axis) are those obtained for the random pack shown in Fig. 6. Experimental results (left vertical axis) were obtained from Chen and Nicholson 11. Although the absolute values of the tortuosities differ, the figure shows that the functional dependence of λ on α are very similar.

Finally, the pore-size distributions of both the regular cubic network (Fig. 7) and random network (Fig. 8) exhibit dramatic changes as a function of the porosity, indicating that there is a nonuniform shrinkage of the vesicle shapes as the porosity is decreased. This is consistent with the behavior proposed by Chen and Nicholson 11 to explain the experimentally observed dependence of the tortuosity on pore volume fraction. In addition, the two-dimensional cross section of the random network shown in Fig. 11 confirms the existence of “lakes” connected, for the most part, by narrow channels.

First-principles modeling studies of the type described in this article provide valuable insight into the geometric properties of the brain's ECS. The networks we considered are constructed by minimizing the configurational energy of the constituent vesicles. We used a simple model free energy, Eq. (1), which is known to provide a good description of the configurational energy of fluid membranes. This is certainly an oversimplification, since it ignores elastic contributions from the cytoskeleton, structural heterogeneities, etc. Nevertheless, we believe that the general features we observe, such as increased tortuosity due to the comparatively large size of the pore void space at network junctions, will remain valid when more sophisticated models for the cells’ mechanical and elastic properties become available. Indeed, our results lend credence to the suggestion by Chen and Nicholson 11 that localized residual space at network junctions, the “lakes,” can help explain both the large tortuosities observed in experiment as well as the characteristic behavior of the tortuosity under osmotic stress. Our results for the random vesicle packs are clearly resolution-limited, and higher resolution studies are required to determine the extent to which models of this type can explain the observed behavior, and resolve remaining questions regarding the importance of viscous or constitutive contributions to the measured tortuosity.