Measuring diffusion using the differential form of Fick's law and magnetic resonance imaging

Diffusion is an important process in many biological and industrial processes. Diffusion coefficients are traditionally measured using integrated solutions of Fick's law for systems with well-defined boundary conditions. We report a simple method for measuring diffusion coefficients in processes without well-defined boundary conditions or without a simple integrated solution using the differential form of Fick's law. Magnetic resonance imaging (MRI) was used to obtain spatially and time-resolved profiles of the diffusion of H2O from an agarose gel to a neighboring D2O reservoir. The differential form of Fick's second law was used to solve for the diffusion coefficient, D=1.3×10−9 m2 s−1. MRI is well suited to this type of analysis as it naturally generates time- and space-resolved images. This analytical method allows for the determination of diffusion coefficients in systems that lack an integral solution to the diffusion equation.


Introduction
Diffusion is a process common to many important porous media systems. The measurement of diffusion coefficients (diffusivity) is central to understanding diffusion processes that occur naturally in porous biological and chemical systems. Numerous techniques and instruments exist for the measurement of diffusivity. The diaphragm cell, the infinite couple and Taylor dispersion are some of the commonly used methods for measuring diffusion coefficients [1]. Nuclear magnetic resonance (NMR) offers benefits for measuring coefficients of self-diffusion, the movement of a tagged solute particle within a chemically identical solvent [1]. The NMR pulsed field gradient (PFG) method is non-invasive and no concentration difference between the two solutions is necessary for the determination of a diffusion coefficient [1].
Diffusive mass transfer in porous media is frequently measured in experiments where one controls or imposes well-known boundary conditions on samples and processes of interest. Ideal boundary conditions permit the integration of Fick's law and fitting of either time or spatially varying concentration data to determine diffusion coefficients [2]- [5]. However, many practical mass transfer processes do not correspond to convenient boundary conditions and therefore have no integrated solution to Fick's law. In other cases, an integrated solution to Fick's law may exist but may be complex and challenging to fit in order to solve for the diffusion coefficient. In either of these situations an alternative approach to measuring the diffusion coefficient is through the differential form of Fick's second law.
The usefulness of the differential approach has been demonstrated by Karger et al [6] in their recent study of the distribution of guest molecules in nanoporous host materials. Interference microscopy was used to generate concentration profiles and determine diffusivity, but the method was limited to optically transparent samples. Magnetic resonance imaging (MRI) is a non-invasive method that can be used to study both transparent and opaque samples. In this work, we use MRI to generate spatially and time-resolved profiles in which the measured signal is proportional to the concentration of H 2 O. This naturally allows for easy use of the differential diffusion equation to determine the diffusion coefficient, D.

Theoretical background
Diffusion can be described by Fick's second law of diffusion in terms of concentration, C, time, t, and space, x, where D is the diffusion coefficient [2]. Knowing ∂C/∂t and ∂ 2 C/∂ x 2 , one can easily solve for D. In order to determine reliable values, it is necessary to smooth the concentration data [7,8]. MRI is a powerful analytical technique but suffers from low intrinsic sensitivity. Without smoothing, the derivatives are unreliable due to a poor signal-to-noise ratio (SNR). The ∂C/∂t and ∂ 2 C/∂ x 2 values can be determined through differentiation with weighted differences and weighted differences of differences equations, respectively, as shown in (2)-(4) [8]. These equations are normalized incremental versions of Gasser's optimal kernel (1,3) and (2,4) functions for continuous functions. The equations do not rely on additional smoothing because the smoothing is implicit in the method. A sample derivation of (2)-(4) is provided in appendix A. The time difference between two consecutive profiles is represented by t. In this work, t = 360 s. The central concentration value is represented by C i , and x is the distance between two adjacent points on a profile. In this work, x = 6.59 × 10 −4 m.
In this work, the concentration values, C, are taken from the MRI profiles, where the measured signal is proportional to the H 2 O concentration. The time and spatial values are naturally generated through MRI. The MRI method used in this work was single point ramped imaging with T 1 enhancement (SPRITE) [9]. SPRITE was chosen because it is a quantitative density imaging technique. The double half-k-space (DHK) method is one of a collection of centricscan SPRITE techniques that sample the origin of k-space first such that no T 1 weighting is introduced at the k-space origin [10]. DHK SPRITE, figure 1, samples one half of k-space from 0 to +k, waits for five times the relaxation time T 1 and then samples the other half of k-space from 0 to -k. The local image signal intensity, S, is given by  where ρ is the local proton density, t p is the experimental encoding time, T * 2 is the effective spin-spin relaxation time and α is the flip angle of the radio frequency (RF) pulse. It can be seen from (5) that when t p T * 2 and when α is constant over the field of view of the image, the signal is directly proportional to the concentration of 1 H, which in this work is equal to the concentration of H 2 O.

Experimental
An agarose gel, 2%, was formed from a mixture of distilled water and agar powder. The solution was heated to boiling and then poured into a glass cylinder (6.8 cm in length, 3 cm in diameter), which was sealed at both ends with Teflon plugs (0.5 cm in length, 3 cm in diameter), as illustrated in figure 2. The gel was cooled passively to room temperature and then placed in a 2 mM solution of CuSO 4 for 48 h. The gel was removed from the CuSO 4 solution, cut to 3 cm with a razor blade and returned to the glass cylinder. A solution of 2 mM CuSO 4 in 99.9% D 2 O was placed on top of the agarose gel. The sample remained horizontal for the duration of the measurements.

5
The one-dimensional (1D) profiles were obtained using a permanent magnet (8.3 MHz) with a 9.5 cm pole gap and a 20 cm long, 16 turn solenoid RF probe (4.8 cm i.d.) controlled from a Tecmag (Houston, TX) Apollo console. The gradient set was homebuilt [11] and water cooled, driven by Techron (Elkhart, IN) 7780 amplifiers.
The repetition time, TR, was set to 2 ms, and the encoding time, t p , was set to 70 µs. The dwell time was set to 8 µs, and the gradient increment G was set to 0.39 G cm −1 . A total of 128 points were acquired and the field of view was 8.7 cm. The 90 • pulse was 20.5 µs, and a pulsewidth of 3 µs (α = 13 • ) was used for the profile acquisition. A total of 32 averages were collected in 1 min.
Bulk T 1 measurements were carried out with an inversion recovery method. A total of 32 τ values were used, biased toward the earlier delays, and four scans were acquired with a delay of 1 s between scans. Bulk T 2 measurements were carried out with a Carr-Purcell-Meiboom-Gill (CPMG) sequence. A total of 16 scans were acquired with 512 echoes and a 200 µs echo time with a 1 s delay between scans. Bulk T * 2 values were determined by fitting to the free induction decay.
1D profiles were collected at even time intervals ( t) after the addition of D 2 O to the agarose gel at 19 • C. t = 180 s from t = 0 to 100 min, t = 360 s from t = 100 to 292 min and t = 660 s from t = 292 to 511 min.

Results and discussion
1D SPRITE DHK profiles were obtained from 0 to 511 min after the D 2 O reservoir was placed adjacent to the agarose gel in the glass tube (figure 3). The signal intensity in the profiles is directly proportional to the H 2 O concentration. By obtaining time-resolved profiles, we can observe the diffusion of H 2 O from the agarose gel to the D 2 O reservoir. It was assumed that the H 2 O concentration is uniform across the vessel with no wall effects. Relaxation times will change as H 2 O forms HDO when it comes in contact with D 2 O. A simple set of measurements was carried out to demonstrate this fact. A pure H 2 O solution had a measured T 1 = 2.0 s and T 2 = 1.9 s, while a 90% D 2 O/10% H 2 O solution had T 1 = 4.8 s and T 2 = 4.1 s at 15 MHz. The DHK SPRITE measurement is immune to such changes in T 1 and T 2 and therefore the profiles remain unaffected.
The smoothing and differentiation operations, (2)-(4), were used to obtain values for ∂ 2 C/∂ x 2 and ∂C/∂t. In the spatial direction, points along the profiles that experienced a significant change in concentration with position were analyzed. These points had concentrations that were approximately 5% different from their neighboring points. Points 71-77 (approximately 5.1-5.6 cm in figure 3) were used in the analysis. Other points along the spatial direction were not considered. Data from early times were discarded because their profiles were affected by geometric blurring due to the nature of the imaging experiment and not the diffusion of H 2 O. Data from late times were discarded because the selected t did not allow for a significant change in concentration to occur before the next measurement. Diffusion coefficients were determined for each of the seven spatial points at times t = 127 to 158 min. Two weighted differences equations were used on the concentration versus time data to obtain two different sets of diffusion coefficients. Equation (3) used five terms that allowed for more smoothing than (2), which used only three terms. The resulting diffusion coefficients from using (2) to determine ∂C/∂t are shown in the second column of table 1. Equation (3) provided similar but more accurate results due to more smoothing. The resulting diffusion coefficients  (2) and (4). b Evaluated using equations (3) and (4). are shown in the third column of table 1. Equation (4) was used consistently for the ∂ 2 C/∂ x 2 value. The data in table 1 show an increasing trend in the diffusion coefficient with time. This may be due to the inherent blurring of the imaging technique. However, the measured diffusion coefficients are still within the experimental error of the measurement, which is estimated to be 12% based on estimated uncertainties in the spatial and temporal directions. A PFG spin echo NMR measurement [12] was employed to determine self-diffusion coefficients of H 2 O at 19 • C. The temperature was maintained at 19 • C due to the water-cooled The results from the differential form of Fick's law were compared with an integrated solution, shown in (6)- (8), where C 0 is the initial concentration of H 2 O in the gel, L 1 is the length of the D 2 O reservoir and L 2 is the length of the H 2 O gel. The detailed solution is provided in appendix B. It is clearly a less desirable, complicated diffusion case. The solution is based on 1D unsteady diffusion from a gel to a well-mixed reservoir. The D 2 O reservoir was well mixed due to passive convection and sample vibration induced by the vibration of the gradients.
It is evident from (6)-(8) that the integrated solution for the current diffusion problem is complicated and fitting to the equation would be challenging. The raw data from profiles within the time range for optimal diffusion were plotted along with (6) with a D value that provided a close match of the experimental curves. An example is illustrated in figure 4. A diffusion coefficient value of ∼ 2.5 × 10 −9 m 2 s −1 allowed for curves with the most overlap with the experimental data; however, this value is too high based on the self-diffusion data. The differential form of Fick's law provides a more reliable answer. If the diffusion problem is not a 1D problem, then one could apply the same methodologies but with slice selection of a 1D portion of the sample. Alternately, one could simply undertake the imaging in 2D or 3D. MRI is well suited to such analyses because it naturally generates spatially and time-resolved images.

The equation obtained from integrating Gasser's kernel functions is
Apply to y = x : Hence, multiply the original formula by 16 807 16 800 :

Appendix B. Integral solution to Fick's second law for diffusion from a gel to a well-mixed reservoir
Diffusion takes place from a gel to a well-mixed reservoir. The system is 1D, but unsteady. The gel initially contains solute at some concentration C 0 , but the solute diffuses out into the reservoir, increasing the concentration in the reservoir from zero. The concentration is continuous at the interface between the gel and the reservoir. At any given time, the overall mass balance gives the concentration in the reservoir in terms of the concentrations across the gel [16,17]: The concentration in the gel is therefore defined by a partial differential equation (PDE) and a corresponding initial condition and two boundary conditions: Unfortunately, integral boundary conditions are not tractable in separation of variables. It is necessary to define an integral parameter here, in order to proceed: The inverse relationship between the cumulative mass and concentration is simply which allows the original PDE to be converted from concentration into cumulative mass, by first taking the derivative in terms of z

The initial condition and two boundary conditions then become
Separation of variables is simplest if the cumulative mass is broken into a steady solution with the above boundary conditions (the last two conditions), and an unsteady solution with an initial condition and two homogeneous boundary conditions: The PDEs for each part of the cumulative mass are The boundary condition that was initially zero remains zero in both cases, while the derivative boundary condition becomes The solution for the steady cumulative mass is a linear profile, where the coefficient b is eliminated by applying one boundary condition, while a is determined from the other boundary condition, The steady solution also determines the initial condition for the unsteady part Separation of variables must be used to solve for the unsteady part. The cumulative mass breaks into the product of two functions as and the two boundary conditions become The PDE then becomes so that the solutions for the two functions are where the time-varying function has been normalized to unity at zero time, in order to match the initial condition simply later on. The boundary conditions eliminate one constant, and the other boundary condition gives an implicit equation for the eigenvalue λ: This equation is most simply solved in dimensionless form, in terms of a parameter β = λL 1 : The solutions for the dimensionless parameter β give an infinite set of solutions, which can be computed numerically for any given geometry. For example, the solutions for the two cases L 1 = L 2 and L 1 = 0.5L 2 are tabulated below up to the fourth smallest value: Once the dimensionless parameter is known it is easily substituted back into the solutions for the two functions, as λ n = β n L 1 .
In these terms, the overall solution is therefore where the Fourier coefficients can be found from the initial condition as where it is again convenient to work in terms of a dimensionless parameter B n = 2L 1 sin β n L 1 + L 2 cos 2 β n β n .
This second parameter can also be tabulated, again using the two example cases L 1 = L 2 and L 1 = 0.5L 2 up to the fourth smallest eigenvalue: Concentration in the reservoir follows from the cumulative mass at the interface as B n sin β n e −β 2 n (Dt/L 2 1 ) β nBn sin β n e −β 2 n (Dt/L 2 1 ) .
The concentration profiles are plotted on the following pages for the two example cases, but it is not possible to construct a general plot here, due to the dependence on the ratio of the gel length to the reservoir length. However, as is clear from the above example solutions, the higher modes of the solution, which decay fastest, are largely independent of the overall split between the gel and reservoir (note the similarity between the values for the fourth mode, for both parameters). This reflects that the concentration profile near the interface is independent of the influence of the lengths of the gel and reservoir-the higher order modes are effectively those with the shortest spatial influence.