MathematicS In Action

The complex transverse water proton magnetization subject to diffusion-encoding magnetic field gradient pulses can be modeled by the Bloch-Torrey partial differential equation (PDE). The associated diffusion MRI signal is the spatial integral of the solution of the Bloch-Torrey PDE. In addition to the signal, the time-dependent apparent diffusion coefficient (ADC) can be obtained from the solution of another partial differential equation, called the HADC model, which was obtained using homogenization techniques. In this paper, we analyze the Bloch-Torrey PDE and the HADC model in the context of geometrical deformations starting from a canonical configuration. To be more concrete, we focused on two analytically defined deformations: bending and twisting. We derived asymptotic models of the diffusion MRI signal and the ADC where the asymptotic parameter indicates the extent of the geometrical deformation. We compute numerically the first three terms of the asymptotic models and illustrate the effects of the deformations by comparing the diffusion MRI signal and the ADC from the canonical configuration with those of the deformed configuration. Thepurpose of this work is to relate the diffusion MRI signal more directly with tissue geometrical parameters.


Introduction
Diffusion magnetic resonance imaging (diffusion MRI) is an imaging modality that can be used to probe the tissue micro-structure by encoding the incoherent motion of water molecules with magnetic field gradient pulses.The motion during the diffusion-encoding time causes a signal attenuation from which the apparent diffusion coefficient (ADC), and possibly higher order diffusion terms, can be calculated [4,11,22].For unrestricted diffusion, the root of the mean squared displacement of molecules is given by x = √ 2 dim σt, where dim is the spatial dimension, σ is the intrinsic diffusion coefficient, and t is the diffusion time.In biological tissue, the diffusion is usually hindered or restricted (for example, by cell membranes).This deviation from unrestricted diffusion can be used to infer information about the tissue micro-structure.
Using diffusion MRI to get tissue micro-structure information in the mammalian brain has been the focus of much experimental and modeling work in recent years [1,3,6,16,17,18,23,24].The predominant approach up to now has been adding the diffusion MRI signal from simple geometrical components and extracting model parameters of interest.Numerous biophysical models subdivide the tissue into compartments described by spheres, ellipsoids, cylinders, and the extra-cellular space (ECS) [1,3,6,10,12,18,19,23].Some model parameters of interest include axon diameter and orientation, neurite density, dendrite structure, the volume fraction and size distribution of cylinder and sphere components and the effective diffusion coefficient or tensor of the ECS.
There is a gold-standard reference model of the diffusion MRI signal, it is the Bloch-Torrey partial differential equation (PDE) that describes the time evolution of the complex transverse water proton magnetization subject to diffusion-encoding magnetic field gradient pulses.The spatial integral of the solution of the Bloch-Torrey PDE provides a reference model for the diffusion MRI signal arising from the geometry of interest.Because of the high computational cost of solving the Bloch-Torrey PDE in complicated cell geometries, this gold standard model has been used primarily as a "forward model" or "simulation framework", in which one changes the inputs parameters such as cell geometry, intrinsic diffusion coefficient, membrane permeability, and study the resulting changes to the MRI signal.This is in contrast to "inverse models", which are used to estimate the model parameters of interest from the MRI signal, the idea being that the "inverse models" have been formulated in such a way so that the model parameters can be correlated to biological information in the imaging voxel."Inverse models" include the biophysical models cited above.
In [13], we presented SpinDoctor, a MATLAB-based diffusion MRI simulation toolbox that solves the Bloch-Torrey PDE using the Finite Element Method (FEM) and an adaptive time stepping method.In addition to the diffusion MRI signal, the time-dependent apparent diffusion coefficient (ADC) can be obtained from the solution of another partial differential equation, called the HADC model, which was obtained using homogenization techniques.SpinDoctor also provides the numerical solution of the HADC model.SpinDoctor provides a user-friendly interface to easily define cell configurations relevant to the brain white matter.In [9], we presented an add-on module of SpinDoctor called the Neuron Module that enables diffusion MRI simulations for a group of pyramidal neurons and a group of spindle neurons whose morphological descriptions were found in the neuron repository NeuroMorpho.Org [2].
In this paper, we continue the Bloch-Torrey PDE based simulation work to further reveal the relationship between the cellular structure and the diffusion MRI signal in the brain white matter.We analyze the Bloch-Torrey PDE and the HADC model in the context of parameterized deformation mappings, starting from a canonical configuration.The canonical configuration we have in mind is a set of straight parallel axons contained in the extra-cellular space.Our idea is to model realistic axons as spatial deformations of canonical configurations of parallel axons.
To be more concrete, we focus on two analytically defined deformations: bending and twisting.We will derive asymptotic models of the diffusion MRI signal and the ADC where the asymptotic parameter indicates the extent of the geometrical deformation.The purpose of this work is to relate the diffusion MRI signal more directly with tissue geometrical parameters.

Theory
We suppose that one would like to simulate a geometrical configuration of axons enclosed in the extra-cellular space (ECS).Let Ω e be the ECS, Ω in i the ith axon.In this paper, we focus on deformations and do not consider the effects of water exchange between the axons and the ECS.Thus, we analyze each compartment individually.In the case where we talk about a general geometrical compartment, we will use the notation Ω.

Bloch-Torrey PDE
In diffusion MRI, a time-varying magnetic field gradient is applied to the tissue to encode water diffusion.For simplicity, we will assume the interfaces between the axons and the ECS are impermeable, meaning the water exchange between compartments is negligible.Denoting the effective time profile of the diffusion-encoding magnetic field gradient by f (t), and let the vector g contain the amplitude and direction information of the magnetic field gradient, the complex transverse water proton magnetization in a compartment Ω (either an axon or the ECS) in the rotating frame satisfies the Bloch-Torrey PDE: where γ = 2.67513 × 10 rad s −1 T −1 is the gyromagnetic ratio of the water proton, I is the imaginary unit, σ is the intrinsic diffusion coefficient in the compartment Ω.The magnetization is a function of position x and time t, and depends on the diffusion gradient vector g and the time profile f (t).Some commonly used time profiles (diffusion-encoding sequences) are the pulsed-gradient spin echo (PGSE) sequence [22] and the oscillating gradient spin echo (OGSE) sequence [7,8].Here, we will consider the PGSE sequence, with two rectangular pulses of duration δ, separated by a time interval ∆ − δ, for which the profile f (t) is where t 1 is the starting time of the first gradient pulse with t 1 + ∆ > T E /2, and T E is the echo time at which the signal is measured.The Bloch-Torrey PDE needs to be supplemented by interface conditions.In this paper, we consider the impermeable boundary condition: where n is the unit outward pointing normal vector and Γ = ∂Ω is the boundary of the compartment Ω (an axon or ECS).The Bloch-Torrey PDE also needs initial conditions: where ρ is the initial spin density in the compartment Ω.The diffusion MRI signal is measured at echo time t = T E > ∆ + δ for PGSE.This signal is the integral of M (x, T E ) in all the compartments: (2.5) In a diffusion MRI experiment, the pulse sequence (time profile f (t)) is usually fixed, while g is varied in amplitude (and possibly also in direction).When g varies only in amplitude (while staying in the same direction), S is plotted against a quantity called the b-value.The b-value depends on g and f (t) and is defined as (2.6) For PGSE, by replacing (2.2) into (2.6), the b-value is [22]: The reason for these definitions is that in a homogeneous medium, the signal attenuation is e −σb , where σ is the intrinsic diffusion coefficient.An important quantity that can be derived from the diffusion MRI signal is the "Apparent Diffusion Coefficient" (ADC), which gives an indication of the root mean squared distance travelled by water molecules in the gradient direction u g = g/∥g∥, averaged over all starting positions: From experimental data, the ADC is numerically computed by a polynomial fit of log S(b).

HADC model
In a previous work [21], a PDE model for the time-dependent ADC was obtained starting from the Bloch-Torrey PDE, using homogenization techniques.In the case of negligible water exchange between compartments (low permeability), there is no coupling between the compartments, at least to the quadratic order in g, which is the ADC term.The ADC in compartment Ω is given by where F (t) = t 0 f (t) dt is the integral of time profile and is a quantity related to the directional gradient of a function ω that is the solution of the homogeneous diffusion equation with Neumann boundary condition and zero initial condition: n being the outward normal and t ∈ [0, T E ].The above set of equations, (2.9)-(2.11),comprise the HADC model.

Canonical configuration and geometrical deformations
To reveal the relationship between the geometrical structure and the diffusion MRI signal, we propose to describe white matter fibers as a deformation of a canonical configuration of parallel axons.The two basic types of deformations that we implement in this paper are 1) bending, and 2) twisting.Both types of deformations will be described by a single parameter, called α twist and α bend .The geometrical structure of the white matter fibers will be defined by these two deformation parameters.Twisting around the z-axis with a twisting parameter α twist is defined by (2.12) Bending on the x − z plane with a bending parameter α bend is defined by We plot in Figure 2.1 a geometrical configuration of 10 cylindrical axons and the ECS before and after deformation.

Derivation of asymptotic models on the deformation parameter
The main aim of our paper is to construct appropriate models to describe the relationship between the deformation parameters α twist and α bend and the diffusion MRI signal as well as the ADC.We will expand the solutions of the Bloch-Torrey PDE and the HADC model as asymptotic series in the deformation parameters α twist and α bend .This approach is expected to work well in the regime of small deformations.

Formulation of the PDEs on the canonical configuration
First, we transform the Bloch-Torrey PDE and the HADC model posed on the deformed geometry Ω into PDEs that are posed on the canonical geometry C. Let r be the space variable in the deformed (by bending or twisting) configuration, whose domain is Ω.The coordinate transformation, maps the canonical configuration defined on C to the deformed configuration on Ω: Let J be the Jacobian of T : We define the composite function for the Bloch-Torrey PDE to be N (x, t) : C → R, where and for the HADC model to be η(x, t) : C → R, where We recall that M (r, t) and ω(r, t) are solutions on the deformed domain Ω, thus, N (x, t) and η(x, t) are the solutions of the respective transformed PDEs.
It is easy to show that the transformed diffusion tensor is : For the twist deformation (α = α twist ): and the transformed diffusion tensor is For the bend deformation (α = α bend ): and the transformed diffusion tensor is The transformed Bloch-Torrey PDE in C is then: n being the outward normal to C. The transformed HADC model is:

Asymptotic expansion of HADC
We now expand the solution of the HADC model in the deformation parameters and match the terms to get the first three terms of the asymptotic expansion.
(3.21) Similarly, in the case of the twisting transformation, using (3.8), the transformed Laplacian operator is (α = α twist ): Also, we define : (3.24)So the transformed Laplacian operator acts as the first and the second correction operators for the Laplacian: where k ∈ {bend, twist}.Using (3.9) and (3.10) for the bend deformation and (3.7) and (3.8) for the twist deformation, the right hand side of the boundary condition of (3.15) becomes: where k ∈ {bend, twist}, and We note that in L twist , we do not expand the trigonometrical functions and keep α twist in the expression.This is because if we simulate a geometry containing long axons, then αz is not a small quantity, there will be a large error if we expand the trigonometrical functions.
The left hand side of the boundary condition of (3.15) becomes: where k ∈ {bend, twist}.For the bend deformation: and for the twist deformation: Finally, we obtain the following equations after matching the terms α j k , with j = 0, 1, 2 and k ∈ {bend, twist}.
For α 0 k , we get the solution of the HADC on the canonical configuration: For α 1 k , we get a PDE that depends on the solution of the previous equation, η 0 : For α 2 k , we get a PDE that depends on the solutions of both of the above PDEs: By using (3.7) and (3.8), we get : The transformed Laplacian operator ∇ • β∇ here is identical to the case of HADC asymptotic expansion.The Iγf (t)g • T (x) operator becomes: where k ∈ {bend, twist}, and (3.45) For the same reason as indicated in previously, we do not expand trigonometrical functions in α twist in the case of the twist deformation.
The left side of boundary condition of (3.11) is also identical to the (3.29).
For simplicity of notation, we define the Bloch-Torrey operator BT := −∇σ∇ + Iγf (t)g • x.We obtain the following equations after matching for α j k , with j = 0, 1, 2, and k ∈ {bend, twist}: For α 0 k , this is the solution of the Bloch-Torrey PDE on the canonical geometry C: For α 1 k , the solution depends on the solution of the above PDE, N 0 : For α 2 k , the solution depends on the solutions of both of the above PDEs:

Numerical implementation
The numerical computations of the asymptotic expansions are done using the diffusion MRI simulation toolbox SpinDoctor [13].We use SpinDoctor to create the geometries, generate finite element (FE) meshes and compute the orders 0, 1, and 2 asymptotic expansions: Firstly, we use SpinDoctor to create a canonical geometry, containing several straight cylindrical axons parallel to the z-axis and an extracellular space wrapped around the axons.Then a finite element mesh is generated for the canonical geometry.The deformed geometries will have finite element meshes that are the analytical deformations of the canonical finite element mesh, described in (2.12) and (2.13).
Since we assumed that the water exchange is negligible, there is no coupling between any compartments for both the Bloch-Torrey PDE and the HADC, and the PDEs are solved independently in each compartment.The finite element discretization is based on continuous piecewise linear basis functions (the P1 finite elements), with a numerically efficient implementation from [20].The time stepping is done automatically using the MATLAB built-in ODE solver ode15s.
Further details about the finite elements matrices construction are contained in the Appendix.

Numerical results
The numerical validation of the asymptotic expansions of the Bloch-Torrey PDE and the HADC model will be conducted in this section.The geometry we use is composed of 10 cylindrical axons and a tightly wrapped ECS, as depicted in Figure 2.1.The radii of the axons are between 2 µm and 3 µm, the exterior width of the ECS is around 40% of the average axon radius, and the height of all the compartments is 20 µm.The diffusion coefficients are set to σ axon = σ ecs = 2 × 10 −3 mm 2 /s and the permeability coefficient set to κ = 0 m/s.The gradient sequence is PGSE (δ = 5 ms, and ∆ = 10 ms).
The reference values are either the ADC obtained by solving the HADC model (Eq.(2.11)) on the deformed geometry Ω or the diffusion MRI signal obtained by solving the Bloch-Torrey PDE (Eq.(2.1)) on Ω.Both of these reference values are obtained using SpinDoctor.The error of the asymptotic model is the difference between the reference values and either

HADC model
First we show the effects of bending and twisting in multiple gradient directions for the HADC model.Being that η 0 gives the ADC of the canonical configuration, η 1 and η 2 could be considered as two corrections.In all the plots that follow, the ADC is normalized by the intrinsic diffusion coefficient σ.In Figure 4.1 we show η 0 , η 1 and η 2 in multiple gradient directions in 3 dimensions.We can see that η 1 provides maximal correction along the z direction.On the other hand, η 2 provides maximal correction along the x − z plane for the bend deformation, and along the x − y plane for the twist deformation.
For the clarity of display, we show further results, which concern the accuracy of our asymptotic model, using two dimensional plots, where a uniform distribution of gradient directions is taken from the x − z plane (y = 0).The reference value is the ADC obtained by solving the HADC model on the deformed geometry Ω.The error of the asymptotic model is the difference between η 0 + η 1 + η 2 and the reference value.In Figure 4.2, we show four curves: the reference value, the asymptotic model (η 0 + η 1 + η 2 , the second order approximation), η 0 (the ADC from the canonical geometry, the zeroth order approximation), and η 0 + η 1 (the first order approximation).We see that frequently, the first order correction is an overcorrection on η 0 and that our second order correction brings the result closer to the reference value.As the deformation parameter increases, the difference between our asymptotic model and the reference value increases, as expected.We note that even though η 0 is the same function on the canonical geometry for both the bend and twist deformations, when transformed to the deformed geometry, its contribution to the ADC is different depending on the specific deformation.This means the computed zeroth order ADC is different for each deformation despite the fact that η 0 is the same function on the canonical geometry C. In Figure 4.3, we show the relative errors of the 0 th , 1 st and 2 nd order approximations, normalized by the reference values.At α bend = 0.05, the maximum 2 nd order approximation error is 1 percent, and the maximum 0 th order approximation error is 14 percent.At α bend = 0.10, the maximum 2 nd order approximation error is 15 percent, and the maximum 0 th order approximation error is 45 percent.At α twist = 0.05, the maximum 2 nd order approximation error is 2 percent, and the maximum 0 th order approximation error is 10 percent.At α twist = 0.10, the maximum 2 nd order approximation error is 22 percent, and the maximum 0 th order approximation error is 45 percent.
Next, we show in Figure 4.4 the relative errors for the axons compartment and for the ECS separately.In general, the axons compartment is much less accurately modelled than the ECS compartment (which is more isotropic).At α bend = 0.05, the maximum axons error over the gradient directions is 8 percent, the maximum ECS error is 1 percent.At α bend = 0.10, the maximum axons error is 40 percent, the maximum ECS error is 8 percent.At α twist = 0.05, the maximum axons error is 16 percent, the maximum ECS error is 1 percent, and at α twist = 0.10, the maximum axons error is 65 percent, the maximum ECS error is 10 percent.

Bloch-Torrey PDE
Now we validate our asymptotic model for the Bloch-Torrey PDE in the same geometries.In Figure 4.5 we show the real part of the normalized signal at b = 3000 s/mm 2 in the canonical geometry, as well as in the bend and twist deformed geometries.In Figure 4.6, we show the relative errors between the 0 th , the 1 st , the 2 nd order approximations and the reference value, for b = 1000 s/mm 2 .For the bend deformation, the 0 th and the 1 st order approximations are indistinguishable in their real parts.For the twist deformation, the 1 st order approximation actually has a higher error than η 0 .
In Figure 4.7, we show the relative errors of the asymptotic model for the axons compartment and for the ECS separately, for b = 500 s/mm 2 and b = 1000 s/mm 2 .For both α bend = 0.05 and α twist = 0.05, the axons errors are larger than the ECS errors.

Convergence order of the asymptotic models
Finally, we show the convergence order of the asymptotic models.In Figure 4.8, we show the relative errors in the ADC of η 0 , η 0 + η 1 , and η 0 + η 1 + η 2 , as α bend and α twist decrease.We see a convergence order of 3, O(α 3 ), for our asymptotic model.
In Figure 4.9, we show the relative errors in the signal of N 0 , N 0 + N 1 , and N 0 + N 1 + N 2 , as α bend and α twist decrease, for b = 500 s/mm 2 .Again, we see a convergence order of 3, O(α 3 ), for our asymptotic model.

Discussion
In the previous section we have shown the accuracy levels of the second order asymptotic models for four geometrical deformations.From Figure 2.1 we can see that at the two smaller deformation values, α bend = 0.05 and α twist = 0.05, there are already visually significant deformations compared to the canonical geometry.It seems that this range of values is sufficient to model significant deviations from straight cylinders and is therefore biological relevant to describe the geometry of the brain white matter.At the higher values that we simulated, α bend = 0.10 and α twist = 0.10, the asymptotic models resulted in much higher errors, but by visual inspection, this larger range of values seems beyond the level of geometrical deviations from straight cylinders that we can expect in the brain white matter.
We have shown that for biologically relevant geometrical deviations, the ADC and the diffusion MRI signal are accurately described as the sum of a zeroth order value (signal or ADC from the straight cylinders) and two orders of corrections.We showed that a first order correction is not sufficient to improve on the zeroth order model, at least two orders of corrections are needed to significantly improve on the zeroth order model.With the second order corrections, the asymptotic models are third order accurate in the geometrical deformation parameters.In addition, the model errors were shown to come mainly from the axons, with the errors from the ECS compartment a much smaller source of error.
This work uses similar mathematical tools as several previous papers focused on the mathematical analysis of the Bloch-Torrey PDE subject to geometrical deformations.In [5], a new mathematical model of Bloch-Torrey PDE in moving and deforming media was introduced.In [14], a rigorous mathematical formalism was introduced to quantify the effect of macroscopicscale tissue motion and deformation in cardiac diffusion MRI.In [15], a new model of the ADC of cardiac diffusion MRI was formulated in the presence of microscopic-scale tissue motion and deformation.
The purpose of this work is to contribute to relating the diffusion MRI signal more directly with the tissue geometrical parameters, the idea being that the diffusion MRI signal and ADC differences between nearby voxels and regions of interest can be modeled by second order corrections due to geometrical deformations with respect to a canonical configuration of straight white matter fibers.Even though the two correction terms we described in this paper are in the forms of partial differential equations and hence are complicated to solve, an intriguing possible future direction is the use of machine learning algorithms to directly map diffusion MRI signals to some geometrical deformation parameters relevant to the white matter fibers in the regions of interest.

Conclusion
We analyzed the Bloch-Torrey PDE and the HADC model in the context of geometrical deformations starting from a canonical configuration, focusing on two analytically defined deformations, bending and twisting.We derived asymptotic models of the diffusion MRI signal and the ADC where the asymptotic parameter indicates the extent of the geometrical deformation.We computed numerically the first three terms of the asymptotic models, the zeroth order model based on the canonical configuration, and two orders of corrections.We showed that a first order correction is not sufficient to improve on the zeroth order model, at least two orders of corrections are needed to significantly improve on the zeroth order model.With the second order corrections, the asymptotic models are third order accurate in the geometrical deformation parameters.In addition, the model errors were shown to come mainly from the axons, with the errors from the ECS compartment a much smaller source of error.The purpose of this work is to contribute to relating the diffusion MRI signal more directly with the tissue geometrical parameters, the idea being that the diffusion MRI signal and ADC differences between nearby voxels and regions of interest can be modeled by second order corrections due to geometrical deformations with respect to a canonical configuration of straight white matter fibers.

Appendix
In this Section, we give details about the construction of the finite elements matrices needed to implement the second order corrections.
For the asymptotic expansion of HADC model, the FE matrices are generated for each compartment with P1 basis function, which are denoted as φ i , for i = 1, . . ., N v , where N v is the number of mesh nodes.The approximate solution of each order of expansion is Nv ξ j i φ i , where ξ j i is the weight of corresponding P1 basis function and j is the order.The zeroth order term is identical to regular HADC model.The mass, stiffness and flux matrices M, S and Q are defined as follows: Applying Green's Theorem for partial derivatives, the integral of correction terms K k,1 and K k,2 multiplying φ k could be written as an volume integral of C and a surface integral of ∂C.The later items cancel G k,1 and G k,2 .Therefore, for the bend deformation, the two correction matrices are: For the twist deformation, the two correction matrices are: 3) The flux matrix for solving (3.35) on the bend geometry is: where u gz and n z (x) are the projections of u g and n(x) onto z-axis, respectively.The flux matrix for solving (3.35) on the twist geometry is: where u gx and u gy are the projections of u g onto x-axis and y-axis, respectively.For numerical computing the solution η 0 of (3.32), the semi-discretized equation is: where ξ 0 i φ i is the approximation of η 0 , ξ 0 i is the coefficient of base functions and ξ = F (t) • 1. where k ∈ {bend, twist}.System of semi-discretized equations (6.6), (6.7), (6.8) can be assembled into one equations as below: where ξ a = [ξ 0 ; ξ 1 ; ξ 2 ] and 0 is the all-zeros matrix with the same dimension as M. (6.9) is solved by MATLAB built-in ODE solver ode15s.This solver will automatically determine the time stepping.When the computation is finished, we decompose ξ a into ξ 0 , ξ 1 and ξ 2 .The approximation of η all equals to the sum of all components above: where k ∈ {bend, twist}.(6.10) Then, ADC coefficient is computed according to equation 2.9.It is worth mentioning that the approximation of η is computed on canonical geometry, but the integration in (2.9) should be performed over deformed coordinates.
The asymptotic expansion of the Bloch-Torrey PDE could be discretized similarly as for the HADC model, the only difference is the diffusion encoding gradient matrices, which are described as:
The FE mesh nodes will be deformed analytically by a coordinate transformation.The deformation process and mesh generation are realized by SpinDoctor routines.

.40) 3 . 3 .
Asymptotic expansion for Bloch-Torrey PDE in the deformation parametersSimilar to the asymptotic expansion of HADC model, we write the solution N (x, t) of (3.11)-(3.13)as a three term expansion:

Figure 4 . 4 .
Figure 4.4.The relative ADC errors between the reference solution and the asymptotic model in 60 gradient directions in the x − z plane (y = 0), in all compartments (blue line), in the axons (red line), and in the ECS (yellow line).The labelled values on the gray circles are given in percent.The displayed angle (from 0 to 360 degrees) is the angle between positive x-axis and the diffusion gradient direction.The angles in bold styles indicate the encoding gradient directions where ADC errors reach minimum.The gradient sequence is PGSE (δ = 5 ms, and ∆ = 10 ms).The ratio of volume of axons and ECS is around 1:1.5.Top left: α bend = 0.05; Top right: α twist = 0.05.Bottom left: α bend = 0.10; Bottom right: α twist = 0.10.

Figure 4 . 5 .
Figure 4.5.The real part of the normalized diffusion MRI signal at b-value = 3000 s/mm 2 , in 180 gradient-directions, which are uniformly distributed on the sphere.The distances from the origin of the black dots as well as the colors are normalized by S 0 .The gradient sequence is PGSE (δ = 5 ms, and ∆ = 10 ms).The diffusion MRI signals of the canonical configuration (left).The signals of the bend deformation by asymptotic model, with α bend = 0.05 (middle).The signals of the twist deformation by asymptotic model, with α twist = 0.05 (right).

Figure 4 . 6 .
Figure 4.6.The relative signal error between 0 th , 1 st and 2 nd order approximations and reference value in 60 directions gradient-directions in the x − z plane (y = 0).The labelled values on the gray circles are given in percent.The displayed angle (from 0 to 360 degrees) is the angle between positive x-axis and the diffusion gradient direction.The b-value = 1000 s/mm 2 and the gradient sequence is PGSE (δ = 5 ms, and ∆ = 10 ms).The real part of the diffusion MRI signal is normalized by the initial signal S 0 .The blue, red, yellow lines represent η 0 , η 0 + η 1 and η 0 + η 1 + η 2 , respectively.Top left: α bend = 0.05; Top right: α twist = 0.05; Bottom left: α bend = 0.10; Bottom right: α twist = 0.10.

Figure 4 . 7 .
Figure 4.7.The relative signal errors between the reference solution and the asymptotic model in 60 gradient directions in the x − z plane (y = 0), in all compartments (blue line), in the axons (red line), and in the ECS (yellow line).The labelled values on the gray circles are given in percent.The displayed angle (from 0 to 360 degrees) is the angle between positive x-axis and the diffusion gradient direction.The angles in bold styles indicate the encoding gradient directions where signal errors reach minimum.The gradient sequence is PGSE (δ = 5 ms, and ∆ = 10 ms).The ratio of volume of axons and ECS is around 1:1.5.Top: α bend = 0.05.Bottom: α twist = 0.05.Left: b = 500 s/mm 2 .Right: b = 1000 s/mm 2 .

Figure 4 . 8 .
Figure 4.8.The ADC relative error (in percent) vs. deformation angle.The diffusion gradient direction is (1, 0, 0).The black, red and blue markers represent zeroth order, first order and second order approximations, respectively.The lines with the same color are the fitting functions.Left: ADC error vs. bend angle; Right: ADC error vs. twist angle.

Figure 4 . 9 .
Figure 4.9.The signal relative error (in percent) vs. deformation angle.The diffusion gradient direction is (1, 0, 0).The black, red and blue markers represent zeroth order, first order and second order approximations, respectively.The lines with the same color are the fitting functions.Left: Signal error vs. bend angle; Right: Signal error vs. twist angle.

Figure 4 . 10 .
Figure 4.10.The signal relative error vs. b-values.The diffusion gradient direction is (1, 0, 0).The black, red and blue markers represent zeroth order, first order and second order approximations, respectively.The lines with the same color are the fitting functions.Left: Signal error vs. b-value, with α bend = 0.06; Right: Signal error vs. b-value, with α twist = 0.06.