Issue #3-4/2019
O.P.Posnansky
Renormalization Group Method: Enhancement of the Sen-Basser Model of the Brain White Matter
Renormalization Group Method: Enhancement of the Sen-Basser Model of the Brain White Matter
We have developed a renormalization group method in order to explore the influence of a large range of geometrical and physical microparameters on the effective diffusion coefficient. The main element of this method comprises of the renormalization of structural parameters from the micro- to the millimeter scale of the random multi-component white matter of the brain. Our approach enhances the Sen-Basser model of the brain white matter [Sen P. N., Basser P. J., Biophys. J., 89 (2005) 2927] by taking the influence of disorder into consideration and allowing quantitative investigation of the sensitivity of the effective diffusion coefficient to the variations of the dominant set of microparameters.
Теги: effective diffusion coefficien renormalization group method sen-basser model метод ренормализационной группы модель сена-бассера эффективный коэффициент диффузии
INTRODUCTION
Brain white matter can be naturally classified into two domains: the intracellular space (ICS) and the extracellular space (ECS). The ICS consists of a tightly packed composite of axons, glia and their cellular extensions, and the ECS is a jelly-like matrix which separates brain cells. The ECS maintains the essential brain functions, such as intercellular communication, nutrient transport and drug delivery. The ECS represents only ~20% of brain parenchymal volume under normal conditions, but can contract and expand dynamically with changes in the brain cell volume during neuronal stimulation and on administration of osmotic agents. Changes in ECS size and geometry are followed by many pathological conditions indicating ischemia, inflammation, and tumour progression.
Taking into account that the brain is a medium composed of randomly arranged domains of different phases, its effective macroscopical properties can be described by the methods developed for the random heterogeneous materials. Using modern magnetic resonance imaging (MRI) techniques, unique information about the organization of white matter, neuronal pathways and their visualization can be provided. Macroscopic properties can easily be acquired giving an average measure of the underlying microstructure. For example, the ECS volume fraction, p, and the effective diffusion coefficient, Deff, are the two major macroscopic parameters commonly employed to give an average description of brain’s microstructure. Thus the problem can be formulated as follows: using the methods of MRI, it is possible to evaluate the macroscopical diffusion coefficient in biological tissue. Further, the sensitivity of the diffusion coefficient to the volume fraction of the comprising phases, the ratio of the microscopic diffusion coefficients both within the cells and outside, the disordered geometry of domains, and the concentration of protons in different areas of tissue may be evaluated. The understanding of the influence of all these parameters on diffusion and spatial fragmentation is important for clinical applications since their significant difference from average values are hallmarks of pathological brain states giving the possibility to detect disease on its very early stages.
In MRI experiments, diffusion is measured as a part of the attenuating signal due to decoherence of phase shifts of water molecules brought about by translational motion in the presence of the external linear magnetic field. Diffusion causes a fundamental limitation to MRI resolution, but on the other hand can provide a contrast mechanism which makes possible the imaging of the molecular displacement spectrum. It is supposed a perfect gradient coil system and compensated local magnetic inhomogeneities for the unbiased estimation of diffusion parameters. However, in the process of signal decay the relaxation of magnetization characterized by longitudinal, T1, and transversal, T2, times is involved. Relaxation can vary in different domains due to variation in the local magnetisation. Nevertheless, at long diffusion times signal attenuation can be described by a single exponential curve with a specific diffusion coefficient:
S/S0 = exp(–b∙Deff), (1)
where b = γ2G2δ2 (Δ – δ/3) is a factor defined according to the Stejskal-Tanner formula [1] (γ is a gyromagnetic constant, S0 is a constant related to relaxation times, Δ is a diffusion time, and G is a magnetic field gradient characterized by duration δ).
Previous attempts have been made to model the brain’s microstructure. It has been shown that the anisotropic properties of the diffusion coefficient are in good agreement with hindered and restricted models of diffusion, especially in the long time scale regime. These diffusion models are also meaningful for the description of anisotropy and fibre tracking in the white matter of brain tissue. The usual approach taken to investigate the heterogeneous nature of brain tissue is to study fibres of packed cylinders arranged in different types of symmetrical arrays such as a square or a hexagon and use the Maxwell-Garnett effective media approximation [2] or depict tissue as a set of parallel cylinders and use the series-parallel approximation [3]. In such models it is difficult to introduce a highly randomized structure which is very natural in biological materials but is quite the opposite of the case in solid-state crystals. Using the approach to diffusion proposed in this article, it is possible to study disordered biological materials and take into account the influence of the complex geometry and topology on the effective transport (diffusive) properties of the tissue.
Based on the following description of the physical process, it is reasonable to guess that the effective medium theory is applicable. The typical diffusion-weighted EPI (echo planar imaging) MRI sequence allows one to acquire information about the dephasing of water molecule spins with durations ranging from Td = 10 ÷ 60 ms. The self-diffusion coefficient of pure water at 36,6° C is 2,94 ∙ 10–9 m2/s. We believe that the presence of glia cells in the extracellular space is consistent with a hindered diffusion coefficient that is lower than 30% of free water at body temperature and is around 2 ∙ 10–9 m2/s. In addition, intracellular diffusion is strongly influenced by neurofilaments and microtubules and thus less than 75% of free water diffusion, estimated as 0.75 · 10–9 m2/s. Using Einstein’s diffusion equation (<x2> = 2DTd) it is possible to calculate the average translation <x2> of water molecules, which is around 5 ÷ 15 μm. The diameter of the fibres varies by about 3 μm. The average distance between fibres is around 1 ÷ 18 μm. It is easy to understand that the water molecule probes all tissue compartments and the MRI signal gives an averaged or homogenized effect of attenuation. Numerical simulation of the dephasing process is rather challenging and takes a lot computational power. Indeed, if the density of free water is 0.995 g/cm3 at a temperature of 36.6 °C and molar mass is 18.01524 g/mol, then in a voxel (2 mm)3 there are around 270 · 1018 molecules or 4.441 · 10–4 mol. Water molecules of size 3,2 Å experience more than 106 collisions during the travelling distance of 15 μm. The huge number of colliding molecules together with the highly inhomogeneous structure of the environment seriously hampers the direct modelling of the discrete random walk. In these circumstances, the mean diffusion can be estimated either by the conventional effective media schemes [4, 5], supposing superposition principles (Maxwell-Garnett approximation) or dedicated spatial averaging technique [6–9].
In this work, we explore a bridge between local and effective global diffusive transport properties and estimate the sensitivity of the diffusion coefficient to the dominant set of micro-parameters: extracellular volume fraction, extracellular diffusion, myelin-sheath diffusion, axon diffusion, mean size of axon, mean size of myelin-sheath, extra-cellular proton density, myelin-sheath proton density, and axon proton density. The goal is to provide a tool that is useful for the entire composition range with heterogeneities occurring at multiple length scales. In the next section the structural model of random heterogeneous material is explained and illustrated with the help of a square tessellation. This is followed by results on the classical theory of the diffusion and its application to local estimations of diffusivity. Finally, the implementation of the renormalization group (RG) method together with numerical results is discussed in the last section.
STRUCTURAL MODEL OF WHITE MATTER
In this section we explain the basic concepts of the renormalization group method in general terms and apply them to the problem of evaluating the effective diffusion coefficient of brain white matter modelled by a square tessellation filling two-dimensional space.
In electron-microscope images of the histological cross-section of biological tissue (white matter of the brain) it is possible to recognise two distinct sets: myelinated axons and the extracellular, gel-like matrix [10–12]. Defining structural and effective physical (diffusion) parameters of a material with randomly distributed local properties is a rather complicated task. Hence we introduce some simplifying assumptions. Myelinated axons can be treated as cylinders which are homogeneously and isotopically statistically distributed in a fluid basin. The symmetry of the problem allows us to split the system into orthogonal 1d longitudinal and 2d transverse subspaces. Then for the mathematical evaluation of the orthogonal features, the tissue can be mapped on a square tessellation (Fig.1a). We supposed that the black square blocks describe properties of the myelinated axon (ICS) and the white square blocks belong to the fluid matrix (ECS). Thus, the white matter is represented as a bundle of cylinders (myelinated axons) randomly embedded in the matrix and can be treated as a two-phase medium, where one of the phases (fibres) has composite properties.
Anatomical sections of white matter reveal that a myelin sheath of neurons has a finite thickness and a longitudinal structure of axon is built from the following components: myelin sheath, characterized by axonal inner, Ra, and outer, Rm, diameters, and the axonal membrane and nodes of Ranvier (Fig.1b). Diffusion across myelin is hindered by layers of lipid bilayers and separate membrane skin. In the modelling of properties of white matter, it is possible to determine the permeability of neurons as a combined effect influencing all substructures. In Fig.1c the transverse structure of axon is depicted with assigned parameters of ECS diffusion, De, and equilibrium concentration of water molecules, ce; ICS is characterized by axonal, Dа, and myelin sheath, Dm, diffusion coefficients and concentrations cа and cm correspondently.
Consider a two-phase structure with a space-filling regular tessellation, whose individual white and black square blocks represent either ECS or ICS (Fig.2a). Such a tessellation can be treated as a set of Wigner-Seitz primitives. Let the overall volume fraction of ECS (white blocks in Fig.2a) in a composite be p0. Then, a fraction p0 of the white blocks possess an upper bound (U) of diffusivity. ICS occupies fraction (1 – p0) and possess a lower bound (L) of diffusivity. Let D^U and D^L be the local diffusions of pure phases U and L whose properties are defined by linear constitutive equations.
The local probability distribution of the diffusion coefficients in this tessellation is:
, (2)
where δ(D^) is the Dirac function. It is supposed that D^0L = D^L and D^0U = D^U. This tessellation is homogeneous on the scale of an individual white square block, and heterogeneous for the black square block representing axon covered by myelin. Instead of the individual blocks, a representative group of adjacent blocks can be considered to be the primary building block of the tessellation. Such a group, called a renormalization cluster, may contain any conceivable configuration of phases and hence could be homogeneous or heterogeneous depending on whether it contains one phase or both phases. Fig.2b provides possible configurations of a renormalization cluster for the square tessellation. Renormalization clusters can vary in size but are normally chosen to be symmetrically shaped in order to maintain the isotropy of the tessellation at each step in the renormalization process. The cluster shown in Fig.2b is the smallest one suitable for the square tessellation. The number of blocks in a renormalization cluster is Ld, where d is the dimensionality of the tessellation and L is the length scale of the cluster. Since the length scale of the unit block of the tessellation is arbitrary, L is the scaling length in the process. For example, L is equal to 2 for the cluster shown in Fig.2b.
The aim of the renormalization process is to replace the original heterogeneous tessellation by an equivalent homogeneous tessellation whose effective diffusion coefficient is to be determined (Fig.2c). This is achieved by a sequence of averaging operations at the scale of the renormalization cluster so that a new probability distribution is obtained at every step. The morphology is assumed to retain its disordered nature and the probability distribution is always composed of a series of delta functions. After n renormalization steps, the distribution is Pn(D^(r→)) with corresponding tessellations. In general, any tessellation is renormalized by replacing every renormalization cluster by a block whose diffusion coefficient is equivalent to the effective diffusion coefficient of that renormalization cluster. As a result of this transformation, the mean cluster size, ξn+1, is related to ξn by equation:
ξn+1 = ξn / L. (3)
That is, each renormalization step reduces the volume of the tessellation by 1/Ld because Ld block is replaced by a single, equivalent block. If, after n → ∞ steps, ξn→∞ comes arbitrary close to zero, the process is complete. At this stage, the tessellation is homogeneous and takes one of the possible D^Un→∞ or D ^Ln→∞ values, and Pn→∞(D^(r→)) is approximately given by:
Pn→∞(D^(r→)) = δ(D^(r→) – D^eff), (4)
where D^eff = D ^U|n L →∞ is the effective diffusion coefficient of the composite being modelled.
This renormalization procedure, often referred to as decimation, was first proposed for site-percolation by Reynolds et al [13].
In order to carry out the renormalization process, it is necessary to identify every possible arrangement of the black and white phases in the renormalization cluster. The number of possible configurations of x phases in a Ld block renormalization cluster is xLd. As the assignment of phases to block is random, the probability of occurrence of a particular configuration of phases in the renormalization cluster is given by the product of the statistically independent probabilities of occurrence of the phases in the blocks constituting the cluster.
Suppose the geometry of ECS allows one to travel from one side of the renormalization cluster to the opposite one without penetrating the other phase. We apply such a condition on different realizations of two-component system, and we can describe the spanned cluster with a joint conditional finite-dimensional distribution R(p). The probability 1 – R(p) is valid for the complementary case.
The connectedness function R(p) for the "2 × 2" RG cluster can be expressed by the polynomial equation
R(p) = p4 + 4p3 (1 – p) + 2p2 (1 – p)2. (5)
Eq. 5 was constructed according to the following rules:
The probability that four, three and two white square blocks (Fig.2b) are present in spanned cluster is p4, p3 (1 – p) and p2 (1 – p)2, respectively.
Only those configurations are available to represent the connectedness of the opposite sides to characterize the spanned cluster.
The probability of the connectedness is then computed by (Table 1): (a) multiplying the probability that rkn blocks are extant with the number of gk blocks configurations; (b) summing these products over all m block (class) configurations.
The whole procedure is depicted in Fig.2b for the conditional probability function.
After the simplification of the polynomial relation Eq. (5), the general equation becomes:
R(p) = 2p2 – p4. (6)
If we take into account the renormalization properties, then
pn+1 = 2p2n – p4n, (7)
where n is a level of the renormalization cluster. The roots (solutions) or fixed points of the algebraic Eq.(7) that belong to the physically meaningful range (0≤p*≤1) are p* = {0,0.618,1}. The first and third roots represent the pristine and vacant lattices (completely full or empty limits). These limits are referred to as being "stable" or "homogeneous" because the recursive procedure for solving Eq.(7) drives the configuration to those limits through successive reduction of the minority phase. At these points the system exhibits translational symmetry and the "stable" fixed points represent "sinks". The second fixed point p* = {0.618} is labelled as being "unstable" or a "source" of the flux of the renormalization parameter pn.
If the composition is exactly at the unstable fixed point, the mean cluster size is infinitely large. Consequently, an infinitely large number of steps are required to scale the heterogeneous system. Thus, the number of steps required in the renormalization procedure increases with the cluster size as p0 approaches the unstable point.
Connectedness function R(p) and its derivative R'(p) are plotted in Fig.2d. The maximum of R'(p) is approached at the fixed point p* = 0,618. Thin dashed lines connecting the bisector and the connectedness functions are the renormalisation of the volume fraction pn according Eq. (7). Within the vicinity of p* = 0,618, the system yields scaling properties and can be treated as a stochastic fractal. Near the fixed point, the typical linear cluster dimension (or the correlation length) ξ = Ln diverges:
ξ ~ |p – p* |–ν, (8)
and the correlation function exponents may be estimated with the help of Eq. (7)
(9)
if for the process the rule
, (10)
is kept at every level n. The critical exponent estimate ν from Eq. (9) is approximately 1.64. Thus the scaling range of the correlation length of the heterogeneous cluster is within the bounds
b ≤ ξ(p) ≤ l, (11)
where l is a macroscopical sample size. In the region p < p* the correlation length is the typical linear size of a ECS island, confined in surrounding fibre bundles, and when p > p* the correlation length is a typical linear size of grouped fibre bundles embedded in a fluid matrix.
LOCAL DIFFUSIVE PROPERTIES
Because we consider diffusion D^(r→) in the stationary regime with no sources of water molecules then locally at every point of the heterogeneous medium, the Laplace law (conservative and potential)
∇→ D^(r→) ∇→c(r→) = 0 (12)
and Fick’s law (linear constitutive relation)
J→ (r→) = –D^ (r→) ∇→c(r→) (13)
hold together with the continuity Cauchy conditions at the two phases interface ∂Г1 and ∂Г2 for completeness. The normal to surface component of the mass flux J→n→ (r→) at the different phases is
J→n→ (r→) | ∂Г1= J→n→ (r→) | ∂Г2, (14)
and the concentration, c(r→), is
c(r→) | ∂Г1= c(r→) | ∂Г2. (15)
For the geometry described in Fig.1b,c boundary conditions (Eqs.(13,14,15)) have the form:
(16)
Here we assumed that temperature and chemical activity are uniform throughout the object and the driving force of mass transport is only a function of concentration. In the real case the axons are covered with myelin sheath interrupted by nodes of Ranvier. Both of them play a role in the electrical properties of axon and transmit a neuronal signal. Myelin has a very dense structure, which leads to the small diffusivity (0,3 · 10–9 m2/s) compared to those of extra- or intraaxonal fluid. Diffusion in the ECS is hindered due to the presence of glia cells and other tissues and is lower than the diffusivity of free water. This fact may be explained on a basis of the presence of glia in the ECS is consistent with an isotropic medium. The ECS then can be characterized with a spherical tensor, which has equal eigenvalues:
(17)
Cylindrical fibres of ICS possess transversal anisotropy, which is expressed in tensor form as:
(18)
After the local homogenization procedure Fick’s law is defined as
<J→(r→)> = –D^eff<∇→c(r→)>, (19)
where the volume averaged mass flux and concentration, respectively, are:
(20)
on the scales ~ L. For symmetry with respect to rotation about one axis, say the z-axis, it can be shown that the effective diffusion tensor has two independent components:
. (21)
This is referred to as transverse isotropic symmetry and describes the geometry of cylindrical fibres embedded in the ECS.
The longitudinal component of the diffusion tensor can be calculated as a volume average [4]:
. (22)
In Eq.(22) f is a volume fraction of axon within Wigner-Seitz primitive and effective concentration of water molecules, ceff, is [4]
(23)
The transverse properties are defined as [4]:
, (24)
where
, (25)
and
ε = сeDe, ε = сaDa, ε = сmDm
RENORMALIZATION GROUP METHOD
In biological tissue it is possible to recognise singular fibres, small groups of isolated fibres and large bundles of fibres confining fluid islands. It is difficult to treat such a system with two distinct scales, where diffusive properties vary either rapidly or slowly. Instead, the system exhibits all ranges of characteristic sizes of heterogeneities, ξ, and they are within a certain range (Eq. (11)). Then, by the virtue of the "2 × 2" – RG model, the diversity of diffusivity on all scales can be derived and used to estimate the effective transversal diffusion coefficient.
Let every square block of the tessellation (Fig.2a) represent the effective microscopic diffusion coefficient in the local area. The structured set consists of two types of squares with diffusion coefficients, D111 (white square) and diffusion coefficients, D112 (black square), respectively. The occurrence probability of a white square is p0 and the occurrence probability of a black one is (1 – p0). The distribution of the phases with two distinct properties is described by the generalized probability density function Eq. (2).
Suppose that random biological tissue can be modelled according to the scheme depicted in Fig.2c. If n is an iterative step in the building of the hierarchical lattice, then for every iterative step the probability density distribution function (Eq. (2)) takes the form:
. (26)
The value Def(U f, j) is an averaged, effective diffusion coefficient calculated for the case of the ECS. Such an effective diffusion coefficient is regarded as an upper bound for the diffusion coefficient. For the initial step, j = 1, the effective diffusion coefficient is given by Def(Uf, 1) = D111 and D111 is the free extracellular water diffusion coefficient. The coefficient
Def(L f, j) represents the inverse of the situation above and therefore Def(Lf, 1) = D112 where D112 is diffusion in the axons.
For heterogeneous tissue averaged upper and lower bounds of the diffusion coefficient can be estimated locally with the help of the Sen-Basser model (Eqs. (24–26)). The upper and lower bounds could be represented by the white and black squares in the building of "2 × 2" – RG cluster. The effective diffusion coefficient then can be estimated with the recursive algorithm. This algorithm has been previously used for the solution of the degradation problem in the solid-state and polymer materials with the help of Monte-Carlo methods [14]. The upper bound of the effective diffusion coefficient is evaluated as [15]
, (27)
and the lower bound takes the form:
, (28)
In general, the renormalization group algorithm consists of the following steps:
Initialize the upper and lower bounds of the diffusion coefficients, namely Def(U f, j) and Def(L f, j).
Initialize the probability pj to find the white square on the "2 × 2" – RG cluster.
Evaluate the upper, Def(Uf, j + 1), (Eqs.24, 25) and lower, Def(Lf, j + 1), (Eqs.24,25) bounds of diffusion coefficients of the whole "2 × 2" – RG according to the Sen-Basser (SB) model (Table 1).
Evaluate the probability R(pj) of the connectedness of the white blocks (Eq.7).
Return to step 1 if the norms |Def(Uf, j + 1) – Def(Uf, j)| and |Def(Lf, j + 1) – Def(Lf, j )| between the effective diffusion coefficients of the current (j + 1)th and the previous (j)th iterations are more then the prescribed error.
Stop the procedure if the norms |Def(Uf, j + 1) – Def(Uf, j) and |Def(Lf, j + 1) – Def(Lf, j )| between the effective diffusion coefficients of the current and previous iteration are less then the prescribed error. The resulting coefficient corresponds to the effective diffusion, Def(Uf(L), j + 1) = Deff(p), of the highly randomised structure with p = p1. Thus, the calculated transversal diffusivity is D11eff = Deff(p).
During the building of the structure described in Fig.2b and 2c, the following three parameters were renormalized (Table 1):
. (29)
System (29) is a system of nonlinear equations representing a flux of parameters in a six-dimensional parametric space.
NUMERICAL RESULTS AND DISCUSSION
In order to obtain a real picture of diffusivity in biological tissue we factored in the parameters experimentally measured in [16, 17]. These parameters are collected in Table 2. The initial value for the upper bound was normalized to Def(Uf, 1) / D111 = 1 (where D111 is diffusion in ECS) and the lower bound was normalized to Def(Lf, 1) / D111 at the initial iterative step (j = 1). The value Def(Lf, 1) was calculated according to the Sen-Basser model (Eq. (4)), where composite intracellular diffusion coefficient is denoted as D2 = D112 ).
In Fig.3a the evolution of the renormalized transversal effective diffusion coefficient D11eff is depicted during j-steps of the recursive procedure. There the upper and lower bounds of the effective diffusion are presented together with the mean value. The observed diffusion coefficient fluctuates between bounds. In the small scale, ξ > b, the fluctuations are gigantic, but when ξ >> b and ξ → L the system exhibits statistically homogenised properties. This phenomenon clearly appears in Figs.3a,b and comprises the convergence of the upper and lower bounds from left to right with an increasing scale ξ. It can be seen that there are three distinct parts in the evolution of the parameters: I linear; II nonlinear; and III stable. Regime I is a region of instability (demonstrating scaling symmetry), during the iteration process the recursive value moves toward the attraction area III (region with translational symmetry). Regime II can be treated as a separatrix area.
In Fig.3a the initial fractions of ECS were p1 = 0.25 (solid line), and p1 = 0.35 (dashed line). It is evident that the upper and lower bounds for the diffusion coefficient converge to a certain value (i.e. attractor) which depends on the initial fraction of the extracellular water. As can be observed in Fig.3a, b, the larger value of the effective diffusion coefficient corresponds with the larger value of the initial fraction of the phase. For the case where p1 = 0.61 the number of iterations is j = 11, this builds a lattice consisting of ~811 bonds. Such a number of bonds is enough to obtain a physical description of the highly randomised heterogeneous biological tissue. The procedure of iteration of the parameters was stopped in accordance with the condition , where the acceptable norm ε = 10–3 was chosen.
In the Fig.4a, the dependence of Deff, Dt, eff, Dl, eff versus p are demonstrated. Every point in these curves is a result of the RG process. For comparison, the Sen-Basser results (classic) are shown as well. The fraction of the cells is (1 – p). Every point of the curve corresponds to the attractor point after j recursion steps. The values of the initial parameters are indicated in the Table 2. The right and left asymptotes (p →{0,1}) converge to the upper and lower bounds of the effective media approximation, respectively.
The intermediate asymptote is located between the upper and lower bounds of the effective diffusion coefficient and demonstrates a characteristic nonlinearity. The linearity disappears in the very high permeability coefficient limit, where intermittence dominates over the localization effects (not shown in the figure). The nonlinearity is usual for the many diffusive (transport) processes and often is depicted with a logarithmic ratio. The nonlinearity is extreme in vicinity of the threshold concentration p* = 0,618 for two-dimensional case. In the realistic case volume fraction of the extracellular water varies from 0.1 till 0.4. As can be recognised from the modelling in such a range of the variance of the volume fraction, the effective transverse diffusive coefficient is more sensitive to the permeability change rather then volume fraction of the extracellular water. Besides, the transversal case gives an inverse dependence of the permeability. Larger permeability contributes in lower transversal effective diffusivity which can lead to the point where two compartments cannot be distinguished any more.
For practical importance it is useful to know the sensitivity
S = |∂log(Deff) / ∂log(X)| (34)
of the apparent diffusion coefficient calculated as a trace of effective diffusion tensor in Eq. (21):
Deff = (2D11eff + D33eff ) / 3. (35)
Initial values (X in Eq.(34)) are taken from the literature, where they were experimentally obtained and collected in Table 2 [2].
The results variation of microparameters are presented in Fig.4b. Analysis of sensitivity gives an indication of how they are associated with brain development and disease states. A similar problem is considered in [18] where the direct Brownian random walk Monte-Carlo numerical experiments were used. In both cases the nonlinear behaviour of effective diffusion versus extracellular volume fraction is confirmed. The reason for such pronounced nonlinear behaviour is in the increase in probability of the blocking of water molecules with a decrease of the extracellular volume fraction. Nonetheless, neuronal cell can exhibit degradation effects leading to the enhancing of permeability; such an effect reduces the blocking of water molecules. This is clearly demonstrated in Fig.4b.
The design of diffusion imaging experiments which can validate the proposed theory is reported in [19]. The essential point of such experiments is to vary the time scale as an argument of diffusion associated with probing the spatial scale. At short range strong fluctuation is expected due to the heterogeneity of brain white matter. Such fluctuations decrease on long term scales leading to the estimation of the non-biased effective apparent diffusion coefficient.
CONCLUSIONS
In summary, we have introduced a scale averaging effective media iterative algorithm to describe a diffusion problem in random heterogeneous biological tissue. The diffusive (transport) properties were studied in the broad range of microscopic and composite parameters for biological multi-compartment heterogeneous tissue composed of cellular and sub-cellular domains. In these domains, the diffusing molecules may have different diffusion coefficients and densities (concentrations), namely within cells and within the outer medium. The results of the proposed effective media scale-averaging iterative scheme were used to explore the effects of a large range of micro-structural and compositional parameters on the apparent (effective) diffusion coefficient on the iterative lattice comprised of ~811 bonds. A self-consistent modelling framework for predicting the long-term scale effective diffusion coefficient was presented. This framework revealed the strong dependence of the apparent diffusion coefficient on the ratio of the microscopic diffusion coefficients of the phases, their concentrations, and the permeability of the cells. ■
Brain white matter can be naturally classified into two domains: the intracellular space (ICS) and the extracellular space (ECS). The ICS consists of a tightly packed composite of axons, glia and their cellular extensions, and the ECS is a jelly-like matrix which separates brain cells. The ECS maintains the essential brain functions, such as intercellular communication, nutrient transport and drug delivery. The ECS represents only ~20% of brain parenchymal volume under normal conditions, but can contract and expand dynamically with changes in the brain cell volume during neuronal stimulation and on administration of osmotic agents. Changes in ECS size and geometry are followed by many pathological conditions indicating ischemia, inflammation, and tumour progression.
Taking into account that the brain is a medium composed of randomly arranged domains of different phases, its effective macroscopical properties can be described by the methods developed for the random heterogeneous materials. Using modern magnetic resonance imaging (MRI) techniques, unique information about the organization of white matter, neuronal pathways and their visualization can be provided. Macroscopic properties can easily be acquired giving an average measure of the underlying microstructure. For example, the ECS volume fraction, p, and the effective diffusion coefficient, Deff, are the two major macroscopic parameters commonly employed to give an average description of brain’s microstructure. Thus the problem can be formulated as follows: using the methods of MRI, it is possible to evaluate the macroscopical diffusion coefficient in biological tissue. Further, the sensitivity of the diffusion coefficient to the volume fraction of the comprising phases, the ratio of the microscopic diffusion coefficients both within the cells and outside, the disordered geometry of domains, and the concentration of protons in different areas of tissue may be evaluated. The understanding of the influence of all these parameters on diffusion and spatial fragmentation is important for clinical applications since their significant difference from average values are hallmarks of pathological brain states giving the possibility to detect disease on its very early stages.
In MRI experiments, diffusion is measured as a part of the attenuating signal due to decoherence of phase shifts of water molecules brought about by translational motion in the presence of the external linear magnetic field. Diffusion causes a fundamental limitation to MRI resolution, but on the other hand can provide a contrast mechanism which makes possible the imaging of the molecular displacement spectrum. It is supposed a perfect gradient coil system and compensated local magnetic inhomogeneities for the unbiased estimation of diffusion parameters. However, in the process of signal decay the relaxation of magnetization characterized by longitudinal, T1, and transversal, T2, times is involved. Relaxation can vary in different domains due to variation in the local magnetisation. Nevertheless, at long diffusion times signal attenuation can be described by a single exponential curve with a specific diffusion coefficient:
S/S0 = exp(–b∙Deff), (1)
where b = γ2G2δ2 (Δ – δ/3) is a factor defined according to the Stejskal-Tanner formula [1] (γ is a gyromagnetic constant, S0 is a constant related to relaxation times, Δ is a diffusion time, and G is a magnetic field gradient characterized by duration δ).
Previous attempts have been made to model the brain’s microstructure. It has been shown that the anisotropic properties of the diffusion coefficient are in good agreement with hindered and restricted models of diffusion, especially in the long time scale regime. These diffusion models are also meaningful for the description of anisotropy and fibre tracking in the white matter of brain tissue. The usual approach taken to investigate the heterogeneous nature of brain tissue is to study fibres of packed cylinders arranged in different types of symmetrical arrays such as a square or a hexagon and use the Maxwell-Garnett effective media approximation [2] or depict tissue as a set of parallel cylinders and use the series-parallel approximation [3]. In such models it is difficult to introduce a highly randomized structure which is very natural in biological materials but is quite the opposite of the case in solid-state crystals. Using the approach to diffusion proposed in this article, it is possible to study disordered biological materials and take into account the influence of the complex geometry and topology on the effective transport (diffusive) properties of the tissue.
Based on the following description of the physical process, it is reasonable to guess that the effective medium theory is applicable. The typical diffusion-weighted EPI (echo planar imaging) MRI sequence allows one to acquire information about the dephasing of water molecule spins with durations ranging from Td = 10 ÷ 60 ms. The self-diffusion coefficient of pure water at 36,6° C is 2,94 ∙ 10–9 m2/s. We believe that the presence of glia cells in the extracellular space is consistent with a hindered diffusion coefficient that is lower than 30% of free water at body temperature and is around 2 ∙ 10–9 m2/s. In addition, intracellular diffusion is strongly influenced by neurofilaments and microtubules and thus less than 75% of free water diffusion, estimated as 0.75 · 10–9 m2/s. Using Einstein’s diffusion equation (<x2> = 2DTd) it is possible to calculate the average translation <x2> of water molecules, which is around 5 ÷ 15 μm. The diameter of the fibres varies by about 3 μm. The average distance between fibres is around 1 ÷ 18 μm. It is easy to understand that the water molecule probes all tissue compartments and the MRI signal gives an averaged or homogenized effect of attenuation. Numerical simulation of the dephasing process is rather challenging and takes a lot computational power. Indeed, if the density of free water is 0.995 g/cm3 at a temperature of 36.6 °C and molar mass is 18.01524 g/mol, then in a voxel (2 mm)3 there are around 270 · 1018 molecules or 4.441 · 10–4 mol. Water molecules of size 3,2 Å experience more than 106 collisions during the travelling distance of 15 μm. The huge number of colliding molecules together with the highly inhomogeneous structure of the environment seriously hampers the direct modelling of the discrete random walk. In these circumstances, the mean diffusion can be estimated either by the conventional effective media schemes [4, 5], supposing superposition principles (Maxwell-Garnett approximation) or dedicated spatial averaging technique [6–9].
In this work, we explore a bridge between local and effective global diffusive transport properties and estimate the sensitivity of the diffusion coefficient to the dominant set of micro-parameters: extracellular volume fraction, extracellular diffusion, myelin-sheath diffusion, axon diffusion, mean size of axon, mean size of myelin-sheath, extra-cellular proton density, myelin-sheath proton density, and axon proton density. The goal is to provide a tool that is useful for the entire composition range with heterogeneities occurring at multiple length scales. In the next section the structural model of random heterogeneous material is explained and illustrated with the help of a square tessellation. This is followed by results on the classical theory of the diffusion and its application to local estimations of diffusivity. Finally, the implementation of the renormalization group (RG) method together with numerical results is discussed in the last section.
STRUCTURAL MODEL OF WHITE MATTER
In this section we explain the basic concepts of the renormalization group method in general terms and apply them to the problem of evaluating the effective diffusion coefficient of brain white matter modelled by a square tessellation filling two-dimensional space.
In electron-microscope images of the histological cross-section of biological tissue (white matter of the brain) it is possible to recognise two distinct sets: myelinated axons and the extracellular, gel-like matrix [10–12]. Defining structural and effective physical (diffusion) parameters of a material with randomly distributed local properties is a rather complicated task. Hence we introduce some simplifying assumptions. Myelinated axons can be treated as cylinders which are homogeneously and isotopically statistically distributed in a fluid basin. The symmetry of the problem allows us to split the system into orthogonal 1d longitudinal and 2d transverse subspaces. Then for the mathematical evaluation of the orthogonal features, the tissue can be mapped on a square tessellation (Fig.1a). We supposed that the black square blocks describe properties of the myelinated axon (ICS) and the white square blocks belong to the fluid matrix (ECS). Thus, the white matter is represented as a bundle of cylinders (myelinated axons) randomly embedded in the matrix and can be treated as a two-phase medium, where one of the phases (fibres) has composite properties.
Anatomical sections of white matter reveal that a myelin sheath of neurons has a finite thickness and a longitudinal structure of axon is built from the following components: myelin sheath, characterized by axonal inner, Ra, and outer, Rm, diameters, and the axonal membrane and nodes of Ranvier (Fig.1b). Diffusion across myelin is hindered by layers of lipid bilayers and separate membrane skin. In the modelling of properties of white matter, it is possible to determine the permeability of neurons as a combined effect influencing all substructures. In Fig.1c the transverse structure of axon is depicted with assigned parameters of ECS diffusion, De, and equilibrium concentration of water molecules, ce; ICS is characterized by axonal, Dа, and myelin sheath, Dm, diffusion coefficients and concentrations cа and cm correspondently.
Consider a two-phase structure with a space-filling regular tessellation, whose individual white and black square blocks represent either ECS or ICS (Fig.2a). Such a tessellation can be treated as a set of Wigner-Seitz primitives. Let the overall volume fraction of ECS (white blocks in Fig.2a) in a composite be p0. Then, a fraction p0 of the white blocks possess an upper bound (U) of diffusivity. ICS occupies fraction (1 – p0) and possess a lower bound (L) of diffusivity. Let D^U and D^L be the local diffusions of pure phases U and L whose properties are defined by linear constitutive equations.
The local probability distribution of the diffusion coefficients in this tessellation is:
, (2)
where δ(D^) is the Dirac function. It is supposed that D^0L = D^L and D^0U = D^U. This tessellation is homogeneous on the scale of an individual white square block, and heterogeneous for the black square block representing axon covered by myelin. Instead of the individual blocks, a representative group of adjacent blocks can be considered to be the primary building block of the tessellation. Such a group, called a renormalization cluster, may contain any conceivable configuration of phases and hence could be homogeneous or heterogeneous depending on whether it contains one phase or both phases. Fig.2b provides possible configurations of a renormalization cluster for the square tessellation. Renormalization clusters can vary in size but are normally chosen to be symmetrically shaped in order to maintain the isotropy of the tessellation at each step in the renormalization process. The cluster shown in Fig.2b is the smallest one suitable for the square tessellation. The number of blocks in a renormalization cluster is Ld, where d is the dimensionality of the tessellation and L is the length scale of the cluster. Since the length scale of the unit block of the tessellation is arbitrary, L is the scaling length in the process. For example, L is equal to 2 for the cluster shown in Fig.2b.
The aim of the renormalization process is to replace the original heterogeneous tessellation by an equivalent homogeneous tessellation whose effective diffusion coefficient is to be determined (Fig.2c). This is achieved by a sequence of averaging operations at the scale of the renormalization cluster so that a new probability distribution is obtained at every step. The morphology is assumed to retain its disordered nature and the probability distribution is always composed of a series of delta functions. After n renormalization steps, the distribution is Pn(D^(r→)) with corresponding tessellations. In general, any tessellation is renormalized by replacing every renormalization cluster by a block whose diffusion coefficient is equivalent to the effective diffusion coefficient of that renormalization cluster. As a result of this transformation, the mean cluster size, ξn+1, is related to ξn by equation:
ξn+1 = ξn / L. (3)
That is, each renormalization step reduces the volume of the tessellation by 1/Ld because Ld block is replaced by a single, equivalent block. If, after n → ∞ steps, ξn→∞ comes arbitrary close to zero, the process is complete. At this stage, the tessellation is homogeneous and takes one of the possible D^Un→∞ or D ^Ln→∞ values, and Pn→∞(D^(r→)) is approximately given by:
Pn→∞(D^(r→)) = δ(D^(r→) – D^eff), (4)
where D^eff = D ^U|n L →∞ is the effective diffusion coefficient of the composite being modelled.
This renormalization procedure, often referred to as decimation, was first proposed for site-percolation by Reynolds et al [13].
In order to carry out the renormalization process, it is necessary to identify every possible arrangement of the black and white phases in the renormalization cluster. The number of possible configurations of x phases in a Ld block renormalization cluster is xLd. As the assignment of phases to block is random, the probability of occurrence of a particular configuration of phases in the renormalization cluster is given by the product of the statistically independent probabilities of occurrence of the phases in the blocks constituting the cluster.
Suppose the geometry of ECS allows one to travel from one side of the renormalization cluster to the opposite one without penetrating the other phase. We apply such a condition on different realizations of two-component system, and we can describe the spanned cluster with a joint conditional finite-dimensional distribution R(p). The probability 1 – R(p) is valid for the complementary case.
The connectedness function R(p) for the "2 × 2" RG cluster can be expressed by the polynomial equation
R(p) = p4 + 4p3 (1 – p) + 2p2 (1 – p)2. (5)
Eq. 5 was constructed according to the following rules:
The probability that four, three and two white square blocks (Fig.2b) are present in spanned cluster is p4, p3 (1 – p) and p2 (1 – p)2, respectively.
Only those configurations are available to represent the connectedness of the opposite sides to characterize the spanned cluster.
The probability of the connectedness is then computed by (Table 1): (a) multiplying the probability that rkn blocks are extant with the number of gk blocks configurations; (b) summing these products over all m block (class) configurations.
The whole procedure is depicted in Fig.2b for the conditional probability function.
After the simplification of the polynomial relation Eq. (5), the general equation becomes:
R(p) = 2p2 – p4. (6)
If we take into account the renormalization properties, then
pn+1 = 2p2n – p4n, (7)
where n is a level of the renormalization cluster. The roots (solutions) or fixed points of the algebraic Eq.(7) that belong to the physically meaningful range (0≤p*≤1) are p* = {0,0.618,1}. The first and third roots represent the pristine and vacant lattices (completely full or empty limits). These limits are referred to as being "stable" or "homogeneous" because the recursive procedure for solving Eq.(7) drives the configuration to those limits through successive reduction of the minority phase. At these points the system exhibits translational symmetry and the "stable" fixed points represent "sinks". The second fixed point p* = {0.618} is labelled as being "unstable" or a "source" of the flux of the renormalization parameter pn.
If the composition is exactly at the unstable fixed point, the mean cluster size is infinitely large. Consequently, an infinitely large number of steps are required to scale the heterogeneous system. Thus, the number of steps required in the renormalization procedure increases with the cluster size as p0 approaches the unstable point.
Connectedness function R(p) and its derivative R'(p) are plotted in Fig.2d. The maximum of R'(p) is approached at the fixed point p* = 0,618. Thin dashed lines connecting the bisector and the connectedness functions are the renormalisation of the volume fraction pn according Eq. (7). Within the vicinity of p* = 0,618, the system yields scaling properties and can be treated as a stochastic fractal. Near the fixed point, the typical linear cluster dimension (or the correlation length) ξ = Ln diverges:
ξ ~ |p – p* |–ν, (8)
and the correlation function exponents may be estimated with the help of Eq. (7)
(9)
if for the process the rule
, (10)
is kept at every level n. The critical exponent estimate ν from Eq. (9) is approximately 1.64. Thus the scaling range of the correlation length of the heterogeneous cluster is within the bounds
b ≤ ξ(p) ≤ l, (11)
where l is a macroscopical sample size. In the region p < p* the correlation length is the typical linear size of a ECS island, confined in surrounding fibre bundles, and when p > p* the correlation length is a typical linear size of grouped fibre bundles embedded in a fluid matrix.
LOCAL DIFFUSIVE PROPERTIES
Because we consider diffusion D^(r→) in the stationary regime with no sources of water molecules then locally at every point of the heterogeneous medium, the Laplace law (conservative and potential)
∇→ D^(r→) ∇→c(r→) = 0 (12)
and Fick’s law (linear constitutive relation)
J→ (r→) = –D^ (r→) ∇→c(r→) (13)
hold together with the continuity Cauchy conditions at the two phases interface ∂Г1 and ∂Г2 for completeness. The normal to surface component of the mass flux J→n→ (r→) at the different phases is
J→n→ (r→) | ∂Г1= J→n→ (r→) | ∂Г2, (14)
and the concentration, c(r→), is
c(r→) | ∂Г1= c(r→) | ∂Г2. (15)
For the geometry described in Fig.1b,c boundary conditions (Eqs.(13,14,15)) have the form:
(16)
Here we assumed that temperature and chemical activity are uniform throughout the object and the driving force of mass transport is only a function of concentration. In the real case the axons are covered with myelin sheath interrupted by nodes of Ranvier. Both of them play a role in the electrical properties of axon and transmit a neuronal signal. Myelin has a very dense structure, which leads to the small diffusivity (0,3 · 10–9 m2/s) compared to those of extra- or intraaxonal fluid. Diffusion in the ECS is hindered due to the presence of glia cells and other tissues and is lower than the diffusivity of free water. This fact may be explained on a basis of the presence of glia in the ECS is consistent with an isotropic medium. The ECS then can be characterized with a spherical tensor, which has equal eigenvalues:
(17)
Cylindrical fibres of ICS possess transversal anisotropy, which is expressed in tensor form as:
(18)
After the local homogenization procedure Fick’s law is defined as
<J→(r→)> = –D^eff<∇→c(r→)>, (19)
where the volume averaged mass flux and concentration, respectively, are:
(20)
on the scales ~ L. For symmetry with respect to rotation about one axis, say the z-axis, it can be shown that the effective diffusion tensor has two independent components:
. (21)
This is referred to as transverse isotropic symmetry and describes the geometry of cylindrical fibres embedded in the ECS.
The longitudinal component of the diffusion tensor can be calculated as a volume average [4]:
. (22)
In Eq.(22) f is a volume fraction of axon within Wigner-Seitz primitive and effective concentration of water molecules, ceff, is [4]
(23)
The transverse properties are defined as [4]:
, (24)
where
, (25)
and
ε = сeDe, ε = сaDa, ε = сmDm
RENORMALIZATION GROUP METHOD
In biological tissue it is possible to recognise singular fibres, small groups of isolated fibres and large bundles of fibres confining fluid islands. It is difficult to treat such a system with two distinct scales, where diffusive properties vary either rapidly or slowly. Instead, the system exhibits all ranges of characteristic sizes of heterogeneities, ξ, and they are within a certain range (Eq. (11)). Then, by the virtue of the "2 × 2" – RG model, the diversity of diffusivity on all scales can be derived and used to estimate the effective transversal diffusion coefficient.
Let every square block of the tessellation (Fig.2a) represent the effective microscopic diffusion coefficient in the local area. The structured set consists of two types of squares with diffusion coefficients, D111 (white square) and diffusion coefficients, D112 (black square), respectively. The occurrence probability of a white square is p0 and the occurrence probability of a black one is (1 – p0). The distribution of the phases with two distinct properties is described by the generalized probability density function Eq. (2).
Suppose that random biological tissue can be modelled according to the scheme depicted in Fig.2c. If n is an iterative step in the building of the hierarchical lattice, then for every iterative step the probability density distribution function (Eq. (2)) takes the form:
. (26)
The value Def(U f, j) is an averaged, effective diffusion coefficient calculated for the case of the ECS. Such an effective diffusion coefficient is regarded as an upper bound for the diffusion coefficient. For the initial step, j = 1, the effective diffusion coefficient is given by Def(Uf, 1) = D111 and D111 is the free extracellular water diffusion coefficient. The coefficient
Def(L f, j) represents the inverse of the situation above and therefore Def(Lf, 1) = D112 where D112 is diffusion in the axons.
For heterogeneous tissue averaged upper and lower bounds of the diffusion coefficient can be estimated locally with the help of the Sen-Basser model (Eqs. (24–26)). The upper and lower bounds could be represented by the white and black squares in the building of "2 × 2" – RG cluster. The effective diffusion coefficient then can be estimated with the recursive algorithm. This algorithm has been previously used for the solution of the degradation problem in the solid-state and polymer materials with the help of Monte-Carlo methods [14]. The upper bound of the effective diffusion coefficient is evaluated as [15]
, (27)
and the lower bound takes the form:
, (28)
In general, the renormalization group algorithm consists of the following steps:
Initialize the upper and lower bounds of the diffusion coefficients, namely Def(U f, j) and Def(L f, j).
Initialize the probability pj to find the white square on the "2 × 2" – RG cluster.
Evaluate the upper, Def(Uf, j + 1), (Eqs.24, 25) and lower, Def(Lf, j + 1), (Eqs.24,25) bounds of diffusion coefficients of the whole "2 × 2" – RG according to the Sen-Basser (SB) model (Table 1).
Evaluate the probability R(pj) of the connectedness of the white blocks (Eq.7).
Return to step 1 if the norms |Def(Uf, j + 1) – Def(Uf, j)| and |Def(Lf, j + 1) – Def(Lf, j )| between the effective diffusion coefficients of the current (j + 1)th and the previous (j)th iterations are more then the prescribed error.
Stop the procedure if the norms |Def(Uf, j + 1) – Def(Uf, j) and |Def(Lf, j + 1) – Def(Lf, j )| between the effective diffusion coefficients of the current and previous iteration are less then the prescribed error. The resulting coefficient corresponds to the effective diffusion, Def(Uf(L), j + 1) = Deff(p), of the highly randomised structure with p = p1. Thus, the calculated transversal diffusivity is D11eff = Deff(p).
During the building of the structure described in Fig.2b and 2c, the following three parameters were renormalized (Table 1):
. (29)
System (29) is a system of nonlinear equations representing a flux of parameters in a six-dimensional parametric space.
NUMERICAL RESULTS AND DISCUSSION
In order to obtain a real picture of diffusivity in biological tissue we factored in the parameters experimentally measured in [16, 17]. These parameters are collected in Table 2. The initial value for the upper bound was normalized to Def(Uf, 1) / D111 = 1 (where D111 is diffusion in ECS) and the lower bound was normalized to Def(Lf, 1) / D111 at the initial iterative step (j = 1). The value Def(Lf, 1) was calculated according to the Sen-Basser model (Eq. (4)), where composite intracellular diffusion coefficient is denoted as D2 = D112 ).
In Fig.3a the evolution of the renormalized transversal effective diffusion coefficient D11eff is depicted during j-steps of the recursive procedure. There the upper and lower bounds of the effective diffusion are presented together with the mean value. The observed diffusion coefficient fluctuates between bounds. In the small scale, ξ > b, the fluctuations are gigantic, but when ξ >> b and ξ → L the system exhibits statistically homogenised properties. This phenomenon clearly appears in Figs.3a,b and comprises the convergence of the upper and lower bounds from left to right with an increasing scale ξ. It can be seen that there are three distinct parts in the evolution of the parameters: I linear; II nonlinear; and III stable. Regime I is a region of instability (demonstrating scaling symmetry), during the iteration process the recursive value moves toward the attraction area III (region with translational symmetry). Regime II can be treated as a separatrix area.
In Fig.3a the initial fractions of ECS were p1 = 0.25 (solid line), and p1 = 0.35 (dashed line). It is evident that the upper and lower bounds for the diffusion coefficient converge to a certain value (i.e. attractor) which depends on the initial fraction of the extracellular water. As can be observed in Fig.3a, b, the larger value of the effective diffusion coefficient corresponds with the larger value of the initial fraction of the phase. For the case where p1 = 0.61 the number of iterations is j = 11, this builds a lattice consisting of ~811 bonds. Such a number of bonds is enough to obtain a physical description of the highly randomised heterogeneous biological tissue. The procedure of iteration of the parameters was stopped in accordance with the condition , where the acceptable norm ε = 10–3 was chosen.
In the Fig.4a, the dependence of Deff, Dt, eff, Dl, eff versus p are demonstrated. Every point in these curves is a result of the RG process. For comparison, the Sen-Basser results (classic) are shown as well. The fraction of the cells is (1 – p). Every point of the curve corresponds to the attractor point after j recursion steps. The values of the initial parameters are indicated in the Table 2. The right and left asymptotes (p →{0,1}) converge to the upper and lower bounds of the effective media approximation, respectively.
The intermediate asymptote is located between the upper and lower bounds of the effective diffusion coefficient and demonstrates a characteristic nonlinearity. The linearity disappears in the very high permeability coefficient limit, where intermittence dominates over the localization effects (not shown in the figure). The nonlinearity is usual for the many diffusive (transport) processes and often is depicted with a logarithmic ratio. The nonlinearity is extreme in vicinity of the threshold concentration p* = 0,618 for two-dimensional case. In the realistic case volume fraction of the extracellular water varies from 0.1 till 0.4. As can be recognised from the modelling in such a range of the variance of the volume fraction, the effective transverse diffusive coefficient is more sensitive to the permeability change rather then volume fraction of the extracellular water. Besides, the transversal case gives an inverse dependence of the permeability. Larger permeability contributes in lower transversal effective diffusivity which can lead to the point where two compartments cannot be distinguished any more.
For practical importance it is useful to know the sensitivity
S = |∂log(Deff) / ∂log(X)| (34)
of the apparent diffusion coefficient calculated as a trace of effective diffusion tensor in Eq. (21):
Deff = (2D11eff + D33eff ) / 3. (35)
Initial values (X in Eq.(34)) are taken from the literature, where they were experimentally obtained and collected in Table 2 [2].
The results variation of microparameters are presented in Fig.4b. Analysis of sensitivity gives an indication of how they are associated with brain development and disease states. A similar problem is considered in [18] where the direct Brownian random walk Monte-Carlo numerical experiments were used. In both cases the nonlinear behaviour of effective diffusion versus extracellular volume fraction is confirmed. The reason for such pronounced nonlinear behaviour is in the increase in probability of the blocking of water molecules with a decrease of the extracellular volume fraction. Nonetheless, neuronal cell can exhibit degradation effects leading to the enhancing of permeability; such an effect reduces the blocking of water molecules. This is clearly demonstrated in Fig.4b.
The design of diffusion imaging experiments which can validate the proposed theory is reported in [19]. The essential point of such experiments is to vary the time scale as an argument of diffusion associated with probing the spatial scale. At short range strong fluctuation is expected due to the heterogeneity of brain white matter. Such fluctuations decrease on long term scales leading to the estimation of the non-biased effective apparent diffusion coefficient.
CONCLUSIONS
In summary, we have introduced a scale averaging effective media iterative algorithm to describe a diffusion problem in random heterogeneous biological tissue. The diffusive (transport) properties were studied in the broad range of microscopic and composite parameters for biological multi-compartment heterogeneous tissue composed of cellular and sub-cellular domains. In these domains, the diffusing molecules may have different diffusion coefficients and densities (concentrations), namely within cells and within the outer medium. The results of the proposed effective media scale-averaging iterative scheme were used to explore the effects of a large range of micro-structural and compositional parameters on the apparent (effective) diffusion coefficient on the iterative lattice comprised of ~811 bonds. A self-consistent modelling framework for predicting the long-term scale effective diffusion coefficient was presented. This framework revealed the strong dependence of the apparent diffusion coefficient on the ratio of the microscopic diffusion coefficients of the phases, their concentrations, and the permeability of the cells. ■
Readers feedback