# Analysis of multisite 2D relaxation exchange NMR

Maxime van Landeghem, Agnes Haber, Jean-Baptiste d'Espinose de Lacaillerie, Bernhard Blümich

► To cite this version:

Maxime van Landeghem, Agnes Haber, Jean-Baptiste d'Espinose de Lacaillerie, Bernhard Blümich. Analysis of multisite 2D relaxation exchange NMR. Concepts in Magnetic Resonance Part A: Bridging Education and Research, 2010, 36A (3), pp.153-169. 10.1002/cmr.a.20157. hal-03984429

HAL Id: hal-03984429

https://hal.science/hal-03984429

Submitted on 12 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analysis of Multisite 2D Relaxation Exchange NMR

MAXIME VAN LANDEGHEM,<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> AGNES HABER,<fnref n="1" />

JEAN-BAPTISTE D'ESPINOSE DE LACAILLERIE,<fnref n="2" /> BERNHARD BLÜMICH<fnref n="1" />

**ABSTRACT:** Two-dimensional (2D) relaxation exchange nuclear magnetic resonance (NMR) is in many ways similar to 2D frequency exchange NMR, except that the encoding times are comparable to the exchange time. This fact prevents the straightforward analysis of the cross-peak intensities in terms of joint probability densities, and quantitative information and understanding can only be obtained by comparison with simulated spectra. Based on simulations, an explanation is proposed as to why interference between relaxation and exchange may lead to asymmetric 2D exchange maps when exchange occurs between more than two sites. Practically, retro-fitting a simulated data set to an experimental one is shown to allow for the determination of the experimental relaxation and exchange parameters. This point is illustrated by studying a two-site model system consisting of interstitial water exchanging within a pack of spherical silica particles. © 2010 Wiley Periodicals, Inc. Concepts Magn Reson Part A 36A: 153–169, 2010.

**KEY WORDS:** NMR; Laplace inversion; exchange NMR; porous materials

## I. INTRODUCTION

Two-dimensional (2D) relaxation exchange NMR is a form of 2D nuclear magnetic resonance (NMR) where molecular motion and transport are mapped by relaxation exchange. This NMR technique is now gaining momentum in a wide variety of applications involving wetting and drying of liquids in porous

media (1–7). This experiment is analogous to 2D frequency exchange NMR but appears more appropriate for the study of the exchange of liquids between pores because relaxivity is expected to depend on pore size, whereas the chemical shift is not (8).

2D Fourier transform frequency exchange NMR spectroscopy (9) is one of the many forms of correlation NMR (10–12), which reveals dynamic processes on the time scale of the longitudinal relaxation time $T_1$. The principle is simple and elegant. During a first period $t_1$, a chemical group is in a particular chemical environment resulting in a given NMR frequency $\omega_1$. This frequency is encoded by recording the corresponding free induction decay (FID) of the magnetization. Then, it exchanges for another chemical environment of NMR frequency $\omega_2$ during a time lag $t_m$, while the magnetization evolution is frozen by

Received 26 November 2009; revised 8 April 2010; accepted 11 April 2010

Correspondence to: Jean-Baptiste d'Espinose de Lacaillerie; E-mail: jean-baptiste.despinose@espci.fr

Concepts in Magnetic Resonance Part A, Vol. 36A(3) 153–169 (2010)

Published online in Wiley InterScience (www.interscience.wiley.com). DOI 10.1002/cmr.a.20157

© 2010 Wiley Periodicals, Inc.

<footnote n="1">¹Institute of Technical and Macromolecular Chemistry, RWTH Aachen University, D-52056 Aachen, Germany</footnote>

<footnote n="2">²Science et Ingénierie de la Matière Molle, UMR 7615 CNRS UPMC, ESPCI ParisTech, 10 rue Vauquelin, 75005 Paris, France</footnote>

<footnote n="3">³Saint-Gobain Recherche, 39 quai Lucien Lefranc, BP 135, 93303 Aubervilliers, France</footnote>
aligning the magnetization with the direction of the
polarizing field. Finally, the FID is measured during
a period starting at time $t = t_1 + t_m$. Ideally, magnet-
ization is transferred during the mixing time $t_m$ from
one site (with a given magnetic environment identi-
fied by an initial frequency $\omega_1$ measured during the
evolution period $t_1$) to another site (with a different
magnetic environment identified by a final frequency
$\omega_2$, which is measured during the detection period
$t_2$). Such a transfer can arise, for example, from
translational and rotational motion of molecules,
rotations or jumps of chemical groups about bonds,
and from cross relaxation (11). After 2D Fourier
transformation, a correlation map is obtained, which
unravels the transfer pathways. From a series of such
maps acquired with different mixing times, the
exchange rates can be deduced (13).

This type of analysis is popular in NMR spectro-
copy, because many different sites can be discrimi-
nated in terms of their chemical shifts in an NMR
spectrum (14). It remains experimentally demanding,
however, because to measure a chemical shift spec-
trum, a highly homogeneous magnetic field is
needed, which requires a sophisticated magnet. Yet,
the nuclear magnetic response of a molecule to a par-
ticular environment is not only characterized by a
chemical shift but also by a rate of decay. Conceptu-
ally, the chemical shift scale can be replaced by a
relaxation or a diffusion scale, so that the occupancy
of the different sites is obtained not by the structure
of the frequency distribution but instead by the distri-
bution of relaxation rates or diffusion coefficients. In
practice, only three or four different sites can usually
be resolved on a relaxation or diffusion scale, but this
disadvantage is largely compensated by the fact that
magnetization decay measurements can easily be
obtained with simple magnets and inhomogeneous
fields. This concept of relaxation exchange NMR is
not new (15, 16), but it has gained recent popularity
with the development of a stable 2D inverse Laplace
transformation (ILT) algorithm (17), a development
reminiscent of Fourier NMR that gained momentum
only once the fast Fourier Transformation algorithm
had been mastered by the NMR community (5, 18).

Despite the fact that relaxation measurements can
be executed with simple hardware and that they do
not require a homogeneous magnet, 2D Laplace
NMR (1, 19–21) did not develop as early as the spec-
troscopic study of motion by 2D Fourier exchange
NMR. There are two likely reasons for this delay.
First, the inversion of the data from the time domain
into the relaxation rate domain involves the ill-
defined ILT (17, 22), which generates artifacts that
are not as well understood as the ones coming from

the Fourier transformation used to recover the fre-
quency space from the temporal signal in Fourier
NMR. Moreover, the numerical ILT is very consum-
ing in terms of CPU time, thus computing an ILT as
proposed by Venkataramanan et al. (17) has only
recently become possible. Second, in most cases,
spatial exchange and spin evolution cannot be con-
sidered separately in relaxation exchange NMR. This
situation is very different from what is commonly
encountered in Fourier exchange NMR. In the latter
case, the slow exchange limit can often be applied
during the evolution and detection periods, which
allows a straightforward prediction of the cross peak
positions. Furthermore, it is precession that is probed
during the evolution periods, and this precession can
be easily frozen during the mixing time by magnet-
ization storage along the direction of the magnetic
field. Consequently, a direct interpretation of the
cross-peak intensities as conditional, or exchange,
probabilities is often possible in 2D Fourier exchange
NMR. It is this relative simplicity that is responsible
for this technique’s early success.

The interpretation of 2D exchange maps, however,
is much less simple. The reason is that the evolution
times must be long enough to allow relaxation encod-
ing. To keep matters simple, one would ideally want
the exchange rate to be small enough so that
exchange would be negligible during the relaxation
encoding periods. This limit is called the slow
exchange limit. The mixing time, on the other hand,
must be long enough so that exchange occurs signifi-
cantly. The applicability of the slow exchange limit
during the evolution period requires that the system
relaxes faster than it exchanges. This means relaxa-
tion is then necessarily significant during the mixing
period, and thus relaxation and exchange must be
considered simultaneously at least during the mixing
period. Conversely, if the exchange rate was high,
while the mixing period was kept short enough to
avoid relaxation, significant exchange would occur
during the relaxation encoding periods.

In short, the fact that it is relaxation which is
probed during the encoding periods implies that
exchange can never be restricted to the mixing period
and must also be considered during the encoding
periods. Consequently, contrary to the case of 2D
Fourier exchange spectra, 2D relaxation exchange
maps can never be interpreted as pure exchange
maps, that is, exchange probability densities. More-
over, relaxation exchange maps are not always sym-
metric, and often the positions of the peaks are
shifted from the cross coordinates of the apparent
relaxation rates identified on the diagonal (8). Con-
sidering these complications, the aim of this article is
to examine if and how exchange dynamics can be extracted from 2D relaxation exchange maps in the case of multisite exchange.

To address this issue and gain a better understanding about what to expect, 2D relaxation exchange maps were simulated for n-site exchange in the presence of $T_2$ and $T_1$ relaxation during the evolution and mixing periods. Following a review of the theory of the evolution of discrete magnetization components under simultaneous relaxation and exchange, the procedure underlying the generation of 2D exchange data sets is recalled. Then, as a first step, results from the literature (2, 8, 23) are illustrated on simulated 1D $T_2$ data sets to establish that tractable exact analytical solutions can indeed be used to interpret exchange data in the case of two-site exchange. In a second step, it is shown how relaxation and exchange parameters can be extracted from experimental maps of two-site relaxation exchange by comparing experimental and simulated data. As the multisite exchange case is too complex for an analytical treatment, it is shown by simulation that irrespective of ILT artifacts the exchange between three and more relaxation pools may lead to asymmetric 2D exchange maps, whereas two-site 2D exchange maps appear to be always symmetric when square data sets are generated.

## II. THEORY OF n-SITE RELAXATION EXCHANGE

### Time Evolution

The focus of this treatment is on magnetization components as a function of relaxation rates $R = R^{(1,2)} = 1/T_{1,2}$, where the subscripts 1 and 2 refer to the longitudinal and the transverse relaxation rates, respectively. This notation is sufficiently general to apply to $T_1-T_1$, $T_1-T_2$ as well as $T_2-T_2$ exchange NMR.

The NMR signal results from the superposition of all the magnetization components

$$s(t) = \sum_{i} [M_{i}(t) - M_{i}^{\text{eq}}]. \quad [1]$$

In the absence of exchange, the evolution of the magnetization component $M_i(t)$ from the relaxation site i follows the differential equation

$$\frac{d}{dt}[M_i(t) - M_i^{\text{eq}}] = -R_i[M_i(t) - M_i^{\text{eq}}], \quad [2]$$

where $R_i = 1/T_1$ and $M_i^{\text{eq}} = M_i^0$ for longitudinal magnetization and $R_i = 1/T_2$ and $M_i^{\text{eq}} = 0$ for transverse

magnetization. When more than one magnetization component are considered, they are collected in a vector, and Eq. [2] becomes

$$\frac{d}{dt}[\mathbf{M}(t) - \mathbf{M}^{\text{eq}}] = -\mathbf{R}[\mathbf{M}(t) - \mathbf{M}^{\text{eq}}], \quad [3]$$

where $\mathbf{M}$ is the vector of magnetization components and $\mathbf{R}$ is a diagonal matrix with the corresponding relaxation rates. This equation has the familiar solution

$$\mathbf{M}(t) - \mathbf{M}^{\text{eq}} = \exp\{-\mathbf{R}(t - t_0)\}[\mathbf{M}(t_0) - \mathbf{M}^{\text{eq}}]. \quad [4]$$

Then, the total NMR signal is given by

$$s(t) = \sum_{i} [\exp\{-R_{i}(t - t_{0})\}[M_{i}(t_{0}) - M_{i}^{\text{eq}}] + M_{i}^{\text{eq}}]. \quad [5]$$

Considering that the magnetization components are defined by their relaxation rates, if normalized, they actually stand for the relaxation rate distribution of the system under study. It thus appears that the NMR signal is a discrete form of a Laplace integral equation. In theory, the recovery of the distribution of relaxation rates through a formal ILT is thus possible. In practice, this inversion is highly unstable and is not actually performed in the context of NMR. To obtain the distribution of NMR relaxation rates, the solution of the inverse Laplace problem is not obtained through a formal inversion but numerically approximated. Consequently, it must be understood that what is commonly referred to in the NMR literature as ILTs are in fact nonnegative least square (NNLS) fits using regularization functions or truncated singular value decomposition (SVD). For the sake of clarity, this article adheres to this terminology.

If the magnetization components exchange, an off-diagonal exchange matrix, or kinetic matrix $\mathbf{K}$, needs to be introduced in Eq. [3], which mixes the magnetization components

$$\frac{d}{dt}[\mathbf{M}(t) - \mathbf{M}^{\text{eq}}] = -(\mathbf{R} + \mathbf{K})[\mathbf{M}(t) - \mathbf{M}^{\text{eq}}]. \quad [6]$$

Equation [6] still has the familiar solution

$$\mathbf{M}(t) - \mathbf{M}^{\text{eq}} = \exp\{-\mathbf{R}(t - t_0)\}[\mathbf{M}(t_0) - \mathbf{M}^{\text{eq}}]. \quad [7]$$

The exponential operator can be understood as a time evolution operator, which produces the magnetization state at time t given the magnetization state at time $t_0$.
The time evolution operator is evaluated by diagonalizing the matrix sum $\mathbf{R} + \mathbf{K}$ with the help of the rotation matrix $\mathbf{Q}$,

$$ \mathbf{U}(t-t_0) = \mathbf{Q} \exp\{-\mathbf{Q}^{-1}(\mathbf{R}+\mathbf{K})\mathbf{Q}(t-t_0)\}\mathbf{Q}^{-1}. \quad [8] $$

## Content of the Kinetic Matrix

Considering a three-site exchange as an example, the kinetic matrix $\mathbf{K}$ can be written as

$$ \mathbf{K} = \begin{bmatrix} k_{11} & -k_{12} & -k_{13} \\ -k_{21} & k_{22} & -k_{23} \\ -k_{31} & -k_{32} & k_{33} \end{bmatrix}, \quad [9] $$

where $k_{ij}$ is the exchange rate between site $i$ and site $j$. As $\mathbf{K}$ only mixes magnetization components and preserves the total magnetization, detailed mass balance requires that whatever amount of magnetization is lost in one component has to be recovered from other components, i.e.:

$$ \mathbf{K} \mathbf{M}^{\text{eq}} = 0. \quad [10] $$

For example, for magnetization component $M_1$, we can write

$$ k_{11} M_1^{\text{eq}} - k_{12} M_2^{\text{eq}} - k_{13} M_3^{\text{eq}} = 0. \quad [11] $$

Similarly, whatever total amount of magnetization is received by one of the magnetization components from other magnetization components has to be restored to them. For example, for magnetization component $M_1$,

$$ k_{12} M_2^{\text{eq}} + k_{13} M_3^{\text{eq}} = k_{21} M_1^{\text{eq}} + k_{31} M_1^{\text{eq}}. \quad [12] $$

This, with the help of the previous equation, leads to

$$ k_{11} - k_{21} - k_{31} = 0, \quad [13] $$

so that the column sums of $\mathbf{K}$ vanish.

By applying the mass balance to each of the three components, six equations are obtained. Furthermore, the total amount of magnetization exchanged is necessarily conserved. Consequently, of the initial nine parameters $k_{ij}$, only four remain independent. In the general case of n-site exchange, the kinetic matrix is determined by $(n-1)^2$ independent parameters.

Note that the rates $k_{ij}$ apply to the forward kinetics from state $i$ to state $j$, so that the rates $k_{ji}$ apply to the reverse kinetics from state $j$ to state $i$. The associated

equilibrium constant $K_{ij}$, which is defined by the ratio $[M_j]/[M_i]$ of concentrations, is given by

$$ K_{ij} = \frac{[M_i]}{[M_j]} = \frac{k_{ij}}{k_{ji}}. \quad [14] $$

## FID in Response to a Single Excitation Pulse

Given the state of the magnetization at time $t_0$ after an excitation pulse, the transverse magnetization at time $t_1$ later is given by Eq. [7]

$$ \mathbf{M}(t_1 + t_0) = \exp\{-[\mathbf{R}^{(2)} + \mathbf{K}](t_1 + t_0)\}\mathbf{M}(t_0). \quad [15] $$

The apparent formal simplicity of this equation is misleading. As we are dealing with an exponential function of a nondiagonal matrix, the prediction of the decay of the total transverse magnetization is not trivial. It involves the determination of the eigenvalues and eigenvectors of the exponential operator as formally stated in Eq. [8], a task which can be performed analytically only if the number of sites $i$ remains limited and is not possible when dealing with unknown continuous distributions of sites.

## General Analytical Solutions for Two-Site Exchange

For reference, the general solution of Eq. [6] for two-site exchange is recalled (23, 24). As there are only two sites, the magnetization vector can be written explicitly as

$$ \mathbf{M} = \begin{bmatrix} M_a \\ M_b \end{bmatrix}, \quad [16] $$

and the kinetic matrix can be written upon balancing the magnetization as

$$ \mathbf{K} = \begin{bmatrix} k_a & -k_b \\ -k_a & k_b \end{bmatrix}, \quad [17] $$

where $k_b = k_{ab}$ and $k_a = k_{ab} M_b/M_a$. Consequently, the matrix $\mathbf{R} + \mathbf{K}$ can be written as

$$ \mathbf{R} + \mathbf{K} = \begin{bmatrix} R_a^{(1,2)} + k_a & -k_b \\ -k_a & R_b^{(1,2)} + k_b \end{bmatrix}. \quad [18] $$

The eigenvalues $\lambda_+$ and $\lambda_-$ of this matrix are given by (23)
$$ \lambda_{+}^{1,2} = \frac{(R_{a}^{1,2} + R_{b}^{1,2} + k_{a} + k_{b}) + \sqrt{(R_{b}^{1,2} + k_{b} - R_{a}^{1,2} - k_{a})^{2} + 4k_{a}k_{b}}}{2}, \quad [19a] $$

$$ \lambda_{-}^{1,2} = \frac{-(R_{a}^{1,2} + R_{b}^{1,2} + k_{a} + k_{b}) - \sqrt{(R_{b}^{1,2} + k_{b} - R_{a}^{1,2} - k_{a})^{2} + 4k_{a}k_{b}}}{2}. \quad [19b] $$

The rotation matrix **Q** can be deduced from the eigenvectors, allowing the calculation from Eqs. [7] and [8] of the exponential of the **R** + **K** matrix, that is of the evolution operator of the magnetization. Experimentally, the eigenvalues of the **R** + **K** matrix are the coordinates of the peaks of the distribution of the exponential decay rates obtained by the numerical approximate of the solution of the inverse Laplace problem (the so-called “ILT,” as explained in the section “time evolution”). As expected, these coordinates tend toward the relaxation rates when the exchange rates become negligibly small compared to the relaxation rates. This is the slow exchange limit discussed in the preceding section. The intensities of the peaks can be obtained analytically from the rotation matrix **Q** in Eq. [8]. Their complicated expressions need not be reproduced here but can be found in the literature [23]. They correspond to the magnetization values only in the slow exchange limit.

## Evolution Magnetization in $T_2-T_2$ Exchange NMR

As an example, the time evolution of the magnetization is examined for a $T_2-T_2$ exchange experiment with a mixing period during which the magnetization is stored along the direction of the magnetic field and experiences a $T_1$ relaxation (Fig. 1).

The magnetization after the mixing time $t_m$ follows from Eq. [7],

$$ \begin{aligned} \mathbf{M}(t_m + t_1 + t_0) - \mathbf{M}^0 \\ = \exp\{-[\mathbf{R}^{(1)} + \mathbf{K}](t_m)\}(\mathbf{M}(t_1 + t_0) - \mathbf{M}^0). \end{aligned} \quad [20] $$

During the evolution period $t_1$ and the detection period $t_2$, the magnetization is given by

$$ \mathbf{M}(t_1 + t_0) = \exp\{-[\mathbf{R}^{(2)} + \mathbf{K}](t_1 + t_0)\}\mathbf{M}(t_0) \quad [21] $$

and

$$ \begin{aligned} \mathbf{M}(t_2 + t_m + t_1 + t_0) \\ = \exp\{-[\mathbf{R}^{(2)} + \mathbf{K}](t_2)\}\mathbf{M}(t_m + t_1 + t_0). \end{aligned} \quad [22] $$

The last three equations combine into

$$ \begin{aligned} \mathbf{M}(t_2 + t_m + t_1 + t_0) \\ = \exp\{-[\mathbf{R}^{(2)} + \mathbf{K}](t_2)\}[\exp\{-[\mathbf{R}^{(1)} + \mathbf{K}](t_m)\} \\ \times (\exp\{-[\mathbf{R}^{(2)} + \mathbf{K}](t_1 + t_0)\}\mathbf{M}(t_0) - \mathbf{M}^0) + \mathbf{M}^0]. \end{aligned} \quad [23] $$

## Statistical Interpretation of Peak Intensities in the Slow Exchange Limit

In 2D Fourier exchange NMR experiments, the fast exchange limit may apply during the mixing time so that relaxation can be neglected during $t_m$, whereas the slow exchange limit applies during the evolution and detection times. This situation is never encountered in 2D relaxation exchange NMR. Only the following three cases may be encountered. Their corresponding evolution operators are given in Table 1.

1. The exchange is slow during the evolution periods, which means that the time scale of exchange is long compared to the time scale of relaxation, and pure relaxation encoding is possible during the evolution periods. Consequently, there is no fast exchange during the mixing period.

2. The exchange is fast compared to relaxation. Then, the fast exchange limit is satisfied during the mixing time. In this case, the mixing period involves only exchange but not relaxation, whereas the evolution and detection periods encode both, relaxation and exchange.

3. The exchange rate is commensurate with the relaxation rates, and both must be considered during the evolution, mixing, and detection periods.

Only in the first case can the 2D relaxation cross peaks be analyzed in terms of exchange kinetics, although not in a simple manner. Being in the slow exchange limit during evolution means that one can measure the relaxation rate $R_i$ at time $t$ and identify it
<figure><img src="image_2.png" /><figcaption>**Figure 1** Pulse sequence for a T₂–T₂ exchange experiment. Following a π/2 read pulse at time t₀, a CPMG sequence of p₁ π pulses equally spaced by t_E(1) is applied during a first evolution period t₁. The transverse magnetization is encoded according to T₂ during t₁. It is then stored along the direction of the magnetic field with the help of a π/2 pulse where it subsequently relaxes with T₁ during the mixing time t_m. After the mixing time, a π/2 read pulse returns the magnetization to the transverse plane, and its T₂ decay is followed for a second time in the detection period t₂ by recording the amplitude of each echo from a second CPMG sequence of n π pulses equally spaced by t_E(2). This sequence is repeated varying p₂ between 1 and n, resulting in an n² time-domain matrix. After ILT, one obtains a 2D T₂–T₂ correlation map.</figcaption></figure>

with the magnetization of the corresponding site $M_i$. In this case, the slow exchange limit applies to the $t_1$ and $t_2$ encoding periods, and Eq. [23] simplifies to

$$ 
\begin{aligned}
\mathbf{M}(t_2 + t_m + t_1 + t_0) \\
&= \exp\{-[\mathbf{R}^{(2)} t_2]\}[\exp\{-[\mathbf{R}^{(1)} + \mathbf{K}](t_m)\} \\
&\quad \times (\exp\{-\mathbf{R}^{(2)}(t_1 + t_0)\}\mathbf{M}(t_0) - \mathbf{M}^0) + \mathbf{M}^0].
\end{aligned} 
\quad [24] 
$$

As the relaxation matrix $\mathbf{R}^{(2)}$ is diagonal, the magnetization of a particular site $i$ at the end of the evolution period is obtained as

$$ 
\begin{aligned}
M_i(t_1 + t_0) &\sim \exp\{-R_i^{(2)} t_1\} M_i(t_0) \\
&= \exp\{-t_1/T_{2,i}\} M_i(t_0).
\end{aligned} 
\quad [25] 
$$

After the mixing time $t_m$, one obtains from Eq. [7],

$$ 
\begin{aligned}
M_j(t_m + t_1 + t_0) - M_j^0 \\
&= \sum_i \{[\exp\{-[\mathbf{R}^{(1)} + \mathbf{K}]t_m\}]_{ji}(M_i(t_1 + t_0) - M_i^0)\},
\end{aligned} 
\quad [26] 
$$

while similarly to what occurred during the evolution period $t_1$, only $T_2$ is encoded during the detection period $t_2$,

$$ 
\begin{aligned}
M_j(t_2 + t_m + t_1 + t_0) \\
&\sim \exp\{-(R_j^{(2)} t_2) M_j(t_m + t_1 + t_0) \\
&= \exp\{-t_2/T_{2,j}\} M_j(t_m + t_1 + t_0).
\end{aligned} 
\quad [27] 
$$

The resultant time-domain signal is then approximated by

$$ 
\begin{aligned}
M_j(t_2 + t_m + t_1 + t_0) &\sim \exp\{-t_2/T_{2,j}\} \\
&\quad \times \left[ \sum_i \{[\exp\{-[\mathbf{R}^{(1)} + \mathbf{K}]t_m\}]_{ji} \right. \\
&\qquad \left. \times (\exp\{-t_2/T_{2,i}\}M_i(t_0) - M_i^0)\} + M_j^0 \right].
\end{aligned} 
\quad [28] 
$$

It is apparent that after 2D inversion over $t_1$ and $t_2$, cross peaks arise at coordinates $(T_{2i}, T_{2j})$ with amplitudes given by

$$ M_{ij} = [\exp\{-(\mathbf{R}^{(1)} + \mathbf{K})t_m\}]_{ij} M_i(t_0). \quad [29] $$

Therefore, when this slow exchange limit is fulfilled for the encoding and detection periods $t_1$ and $t_2$, the cross-peak amplitude represents the magnetization from site $j$ that originated from site $i$ before the mixing time $t_m$.

For the sake of understanding, the virtual case of fast exchange during the mixing is considered. Then, we would have

$$ \exp\{-(\mathbf{R}^{(1)} + \mathbf{K})t_m\} \sim \exp\{-\mathbf{K}t_m\}. \quad [30] $$

In this limit, magnetization would have been conserved during the exchange and thus, upon normalization, could have been assigned to the probability of site occupancy. The cross peaks would then have mapped the joint probability of molecules or spins being at site $j$ at time $t_0 + t_1 + t_m \sim t_0 + t_m$ and at site $i$ at time $t_0 + t_1 \sim t_0$. This joint probability density can be expressed in terms of the conditional

<table><caption>**Table 1** Evolution Operators of the Three Cases That May Occur in 2D Relaxation Exchange NMR as Described in the Text</caption><thead><tr><th></th><th>U(t<sub>1</sub>)</th><th>U(t<sub>m</sub>)</th><th>U(t<sub>2</sub>)</th></tr></thead><tbody><tr><td>Case (1)</td><td>exp{−[R](t<sub>1</sub>)}</td><td>exp{−[R + K](t<sub>m</sub>)}</td><td>exp{−[R](t<sub>2</sub>)}</td></tr><tr><td>Case (2)</td><td>exp{−[R + K](t<sub>1</sub>)}</td><td>exp{−[K](t<sub>m</sub>)}</td><td>exp{−[R + K](t<sub>2</sub>)}</td></tr><tr><td>Case (3)</td><td>exp{−[R + K](t<sub>1</sub>)}</td><td>exp{−[R + K](t<sub>m</sub>)}</td><td>exp{−[R + K](t<sub>2</sub>)}</td></tr></tbody></table>

The time periods $t_1$, $t_m$, and $t_2$ are defined in Fig. 1. For the T₂–T₂ exchange experiment schematized in Fig. 1, **R** stands for **R**⁽¹⁾ during $t_m$ and for **R**⁽²⁾ during $t_1$ and $t_2$.
probability density of a spin being at site j knowing it started at time i following

$$M_{ji}(t_0 + t_m) = \langle R_j, (t_0 + t_m) | R_i, t_0 \rangle_{ji} M_i(t_0). \quad [31]$$

Consequently, comparing Eqs. [27] and [29], the exponential exchange operator could have been identified with the conditional probability matrix

$$\exp\{-\mathbf{K} t_m\} = [\langle R_j, (t_0 + t_m) | R_i, t_0 \rangle_{ji}], \quad [32]$$

and the cross-peak intensities would have directly mapped the elements of the exponential of the kinetic matrix, or equivalently of the conditional probability, as

$$\begin{aligned} M_{ji}(t_0 + t_m) &= \langle R_j, (t_0 + t_m) | R_i, t_0 \rangle_{ji} M_j(t_0) \\ &= [\exp\{-\mathbf{K} t_m\}]_{ji} M_i(t_0). \end{aligned} \quad [33]$$

This means that the evolution of the intensities of the peaks with the length of the mixing period would have decayed according to the eigenvalues of the exchange matrix.

However, in reality, the fast exchange limit cannot be fulfilled during the mixing time when the slow exchange limit applies during the encoding times. As a result, magnetization is not conserved during the mixing time because of relaxation. Nevertheless, if all the sites had the same relaxation rate $R_{\text{all}}$, one could still write

$$\begin{aligned} &[\exp\{-(\mathbf{R}^1 + \mathbf{K})t_m\}]_{ji} \\ &= \exp(-R_{\text{all}} t_m)[\exp\{-\mathbf{K}t_m\}]_{ji} \\ &= \exp(-t_m/T_{\text{1,all}})[\exp\{-\mathbf{K}t_m\}]_{ji}, \end{aligned} \quad [34]$$

and the cross-peak intensities would still be a representation of the elements of the kinetic matrix albeit weighted by relaxation during the mixing time

$$\begin{aligned} M_{ji}(t_0 + t_m) &= \exp(-R_{\text{all}} t_m) \langle R_j, (t_0 + t_m) | R_i, t_0 \rangle_{ji} M_i(t_0) \\ &= \exp(-R_{\text{all}} t_m) [\exp\{-\mathbf{K} t_m\}]_{ji} M_i(t_0). \end{aligned} \quad [35]$$

Nevertheless, 2D relaxation exchange NMR is of interest only when the sites relax differently, so that the case above is not pertinent. In general, even in the slow exchange limit, one cannot go further other than writing

$$\begin{aligned} M_{ji}(t_0 + t_m) &= [\exp\{-(\mathbf{R}^{(1)} + \mathbf{K})t_m\}]_{ji} M_i(t_0) \\ &= [\mathbf{Q} \exp\{-\mathbf{Q}^{-1}(\mathbf{R}^{(1)} + \mathbf{K})\mathbf{Q} t_m\} \mathbf{Q}^{-1}]_{ji} M_i(t_0). \end{aligned} \quad [36]$$

The other two cases (as explained at the beginning of this section, namely exchange compared to relaxation and exchange commensurate with relaxation) are even more complex as magnetization is not conserved or follows multiple single relaxation paths. The evolution of magnetization must be predicted working back from the general Eqs. [8] and [23]. As a consequence, 2D relaxation exchange maps cannot be interpreted in terms of conditional probability densities, and the observed decay rates are not associated with single-site relaxation rates (8, 23).

## Multisite Relaxation Exchange

In the case of an interconnected system of more than two sites, tracking the general analytical expression of the rotation matrix $\mathbf{Q}$ to diagonalize the evolution operator becomes extremely tedious at best and close to impracticable in most cases. Furthermore, as the slow exchange limit is typically violated in relaxation exchange NMR, 2D exchange maps cannot simply be interpreted in terms of joint probability densities. At this stage of the discussion, as no obvious generalizations can be made from the two-site analytical solution, it becomes clear that only numerical simulations that take relaxation and exchange into account at all time periods of the 2D exchange experiment can support the interpretation of the experimental data.

## III. NUMERICAL SIMULATIONS

The simulation procedure provides the opportunity to generate time-domain exchange data sets free of experimental noise, which can be processed by one of the available so-called ILT algorithms to evaluate methods of extracting relaxation, exchange, and population parameters from the exchange maps. Given a chosen set of $\mathbf{K}$, $\mathbf{R}^{(2)}$, and $\mathbf{R}^{(1)}$ matrices and an initial magnetization vector $\mathbf{M}(t_0)$, which typically is proportional to the population of sites, the magnetization vector $\mathbf{M}(t_1 + t_0)$ is calculated according to Eq. [21] for each site $i$ and at each time $t_1 = p_1 t_E^{(1)}$ with $p_1$ varying between 1 and $n$:

$$\begin{aligned} M_i(t_1 + t_0) &= M_i(p_1 t_E^{(1)} + t_0) \\ &= \sum_{j=1}^{n} ([\exp\{-(\mathbf{R}^{(2)} + \mathbf{K})(p_1 t_E^{(1)})\}]_{ji} M_j(t_0)). \end{aligned} \quad [37]$$

This provides a 1D time-domain data set for each site. Then, for each value of $p_1$ the magnetization vector $\mathbf{M}(t_2 + t_m + p_1 t_E^{(1)} + t_0)$ is calculated according to
<figure><img src="image_3.png" /><figcaption>**Figure 2** Notation for 1D relaxation time distributions of two exchanging sites. [Color figure can be viewed in the online issue, which is available at www.interscience. wiley.com.]</figcaption></figure>

Eq. [23] for each site *j* and at each time *t*₂ = *p*₂ *t*<sub>E</sub><sup>(2)</sup> with *p*₂ running from 1 to *n*. Thus, for each site *j*, an individual *n*² 2D data set is obtained in the time domain,

$$
\begin{align*}
M_j(t_2 + t_m + t_1 + t_0) &= M_j(p_2 t_E^{(2)} + t_m + p_1 t_E^{(1)} + t_0) \\
&= \sum_{r=1}^{n} \{[\exp\{-(\mathbf{R}^{(2)} + \mathbf{K})p_2 t_E^{(2)}\}]_{jr} \sum_{i=1}^{n} ([\exp\{-(\mathbf{R}^{(1)} + \mathbf{K})t_m\}]_{ri} (M_i(p_1 t_E^{(1)} - M_i^0)) + M_r^0)\}. \quad [38]
\end{align*}
$$

Finally, the simulated 2D data set is the sum of the
individual 2D data sets for each site *j*,

$$
\begin{align}
S(t_1, t_m, t_2) &= S(p_2 t_E^{(2)}, t_m, p_1 t_E^{(1)}) \nonumber \\
&= \sum_{j=1}^{n} M_j(p_2 t_E^{(2)}, t_m, p_1 t_E^{(1)}). \tag{39}
\end{align}
$$

This 2D time-domain array simulates the data from
the 2D relaxation exchange NMR experiment (Fig.
3) (16). Finally, the 2D T<sub>2</sub>-T<sub>2</sub> exchange map is
obtained by 2D ILT of this array. The numerical
approximation of the inversion, commonly called
ILT, is performed with a software developed by the
group of Prof. P. Callaghan and based on the method
published by Venkataramanan (namely a NNLS fit
using SVD) (17, 20, 25, 26).

**Procedure for 1D Simulations**

As a preparatory check for 2D exchange data simula-
tions, and to gain prior insight into relaxation
exchange NMR, 1D CPMG data are analyzed. Time-
domain data sets were generated for two-site
exchange using Eq. [39]

$$
\begin{align*}
S_{1D}^{\text{two sites}}(t_1, t_m, t_2) &= S_{1D}^{\text{two sites}}(p_2 t_E^2, t_m, p_1 t_E^1) \\
&= M_A(p_2 t_E^2, t_m, p_1 t_E^1) + M_B(p_2 t_E^2, t_m, p_1 t_E^1) \quad [40]
\end{align*}
$$

and subsequently converted into distributions of
relaxation times by ILT. The notation used in Eq.
[38] is the following (Fig. 2):

* $T_{2\text{app}}^A$ and $T_{2\text{app}}^B$ denote the coordinates of the peaks on a transverse relaxation times scale obtained by ILT of the simulated data.

* $A^A$ and $A^B$ denote the amplitudes of the peaks in the relaxation time distribution of peaks A and B, respectively, calculated from the simulated data by ILT.

* $I_A$ and $I_B$ denote the areas of the ILT peaks of sites A and B, respectively.

The input parameters chosen for the 1D simula-
tions are summarized in Table 2. The exchange time,
e.g., the inverse of the exchange rate, is varied
between 0.01 and 500 ms. For these parameters, 1D
relaxation decays were simulated and inverted by
ILT to extract the T<sub>2app</sub> coordinates of the peaks in
the relaxation time distributions, the peak amplitudes,
and the peak integrals. The results are summarized
graphically in Figs. 3(a,b).

To check for a possible influence of the sampling
interval $t_E$ on the output results, two sampling rates
commensurate with the used exchange time were
tested, namely $t_E = 0.1$ and 0.01 ms. The simulations
led to the same results in both cases. The data pre-
sented here correspond to $t_E = 0.1$ ms.

**Results of 1D Simulations**

The influence of the exchange rate on the peak inte-
grals and peak positions is illustrated in Figs. 3 and
4. When the exchange time becomes short compared
to the relaxation times, one of the peaks grows at the
expense of the other. For an exchange time lower
than 1 ms, only a single peak is observed, and the
peak integral corresponds to the number of spins
from both sites [Fig. 3(b)]. On the other hand, when
the exchange time is at least 10 times longer than the

<table><caption>**Table 2** Input Parameters for Simulations of 1D $T_2$ Distributions</caption>
  <thead>
    <tr>
      <th>Parameter</th>
      <th>Value</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td><i>T</i><sub>2</sub><sup>A</sup></td>
      <td>2.5 ms</td>
    </tr>
    <tr>
      <td><i>T</i><sub>2</sub><sup>B</sup></td>
      <td>25 ms</td>
    </tr>
    <tr>
      <td><i>A</i><sup>A</sup></td>
      <td>0.5</td>
    </tr>
    <tr>
      <td><i>A</i><sup>B</sup></td>
      <td>0.5</td>
    </tr>
  </tbody>
</table>
<figure><img src="image_4.png" /><figcaption>**Figure 3** Results of 1D simulations for two-site exchange. (a) Apparent relaxation times derived from relaxation time distributions and eigenvalues of the **R** + **K** matrix. (b) Peak integrals of relaxation time distributions. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]</figcaption></figure>

longest relaxation time, the peak integrals provide the correct population values. For shorter exchange times, the population of the peak with the longer relaxation time is overestimated and that of the other is underestimated. This observation is in line with the analytical predictions (2, 23).

To verify that these results were not an artifact of the ILT algorithm, the simulated relaxation curves were fitted with two exponential decay functions. The fit results confirmed that the peak integrals of the relaxation time distributions correspond to the amplitudes of the exponential fit functions. The amplitudes follow exactly the same pattern as the peak integrals, validating the use of peak integrals in the analysis of relaxation-time distributions to quantify spin populations.

It is emphasized that the relaxation times extracted from the maxima of the relaxation time distributions are apparent relaxation times, corresponding to the eigenvalues of the sum of the exchange and relaxation matrices and as such are strongly affected by the exchange rate (Fig. 4). An exchange time 100 times longer than the longest $T_2$ component is required to obtain an unbiased estimate of $T_2$. Shorter exchange times lower the apparent transverse relaxation times. When the exchange approaches infinity, the shorter apparent $T_2$ value tends toward zero, whereas the longer one shifts toward a value corresponding to the population-weighted average of the relaxation rates, namely slightly below 5 ms in the present case.

In the intermediate exchange region between the fast and the slow exchange limits, the apparent relaxation times $T_{2\text{app}}$ follows the eigenvalues of the **R** + **K** matrix [Fig. 3(a)]. The excellent agreement between the eigenvalues and the apparent relaxation times shows that, as expected from Eq. [6], the apparent relaxation times are indeed given by the eigenvalues of the **R** + **K** matrix. To anticipate the discussion on 2D exchange experiments, we note that a 1D relaxation time distribution can be understood in terms of a projection of a 2D distribution onto one axis, and thus that the peaks in 2D distributions will be shifted by exchange along the relaxation-time axis in the same way as those in 1D distributions.

The results reported so far have all been simulated with one set of $T_2$ values, i.e., with $T_2^A = 2.5$ ms and $T_2^B = 25$ ms. Similar results were obtained with other choices of relaxation times as input variables.

<figure><img src="image_5.png" /><figcaption>**Figure 4** Examples of 1D spectra of two-site exchange for fast exchange ($k = 10 \text{ ms}^{-1}$), two intermediate exchanges ($k = 0.5$ and $0.1 \text{ ms}^{-1}$), and slow exchange ($k = 0.002 \text{ ms}^{-1}$). [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5** 2D distributions for two-site exchange simulated with the parameters of Table 3 with different exchange rates in the limit of infinite longitudinal relaxation times (a) and with longitudinal relaxation times equal to six times the transverse relaxation times (b).</figcaption></figure>

The dependences of the amplitudes and apparent relaxation times on the exchange time are similar in all cases, which corroborates the conclusions derived above. Moreover, if $T_2$ is replaced by $T_1$, the results of this study do not change. The $T_1$ values are affected by the exchange in the same way as the $T_2$ values.

## 2D Simulations: Two-Site Exchange

The time-domain data for two-site exchange were transformed to 2D exchange maps (Fig. 5) by 2D ILT. The 2D distributions exhibit two diagonal peaks and two off-diagonal peaks, which reveal the exchange between the two sites.

Depending on the exchange time and the exchange rate, one diagonal peak is observed for fast exchange, two diagonal peaks for slow exchange, and additional cross peaks for intermediate exchange. As relaxation during the mixing period complicates the analysis of peak integrals (see Eq. [23]), simulations with the longitudinal relaxation time on ($T_1$ finite) and off ($T_1$ infinite) are compared. Simulations performed in this way with the parameters specified in Table 3 produced the exchange maps of Fig. 5.

The mixing time was set equal to the exchange time in each case except in the first case (top).

When the exchange is very fast (top left), the exchange map exhibits a single peak in accordance with the 1D study. When the exchange time is increased and becomes comparable to the relaxation times, cross peaks appear as well as two distinct diagonal peaks (top right, bottom left). When the exchange time becomes long, i.e., about 10 times the longest value of $T_2$, the exchange is too slow to be seen and the cross peaks disappear (bottom right).

The peak intensities corresponding to the peak integrals in the 2D maps vary with the mixing time $t_m$. The diagonal peaks decay (Fig. 6, top) and the cross peaks initially grow exponentially with the mixing time $t_m$ (Fig. 6, bottom). This has been simulated for two-site exchange with $T_2 = 1$ ms and 10 ms for the two sites (Fig. 6). Without relaxation during the mixing time, the cross-peak intensity grows exponentially toward a dynamic equilibrium value (Fig. 6, bottom left). The time constants of the cross-peak build-up were determined via exponential fits. They are found to correspond to the inverse of the nonzero eigenvalue of the exchange matrix $K$, i.e., to the input parameters $1/k_1 = 50$ ms and $1/k_2 = 100$ ms used in the simulations. As longitudinal relaxation has been turned off in the simulation, any change of peak intensity can indeed only be due to the exchange, so that the peak intensities can only change by mixing between the two populations.

When longitudinal relaxation is included during the mixing time $t_m$, the cross-peak intensity goes to zero for long mixing times (Fig. 6, bottom right), and only when the longitudinal relaxation times are of the order of the exchange time and much larger than the transverse relaxation times did the cross peaks become intense. The build-up of the cross-peak intensity is no longer determined by the nonzero eigenvalue of the exchange matrix $K$, as its growth is now affected by the $T_1$ relaxation of both sites.

<table><caption>**Table 3** Parameters Used to Produce the Exchange Maps of Fig. 7</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>$T_2^A$</td><td>2.5 ms</td></tr><tr><td>$T_2^B$</td><td>25 ms</td></tr><tr><td>$T_1^A$</td><td>Infinite, 15 ms</td></tr><tr><td>$T_1^B$</td><td>Infinite, 150 ms</td></tr><tr><td>$k_1$</td><td>10 ms<sup>-1</sup></td></tr><tr><td>$k_2$</td><td>0.1 ms<sup>-1</sup></td></tr><tr><td>$k_3$</td><td>0.01 ms<sup>-1</sup></td></tr><tr><td>$k_4$</td><td>0 ms<sup>-1</sup></td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Figure 6** Evolution of peak integrals with the mixing time $t_m$ in 2D maps for two-site exchange. Top: Integrals of the diagonal peaks with the short (left) and the long (right) relaxation time with $T_1$ relaxation during the mixing time. Bottom: integral of cross peaks without (left) and with (left) longitudinal relaxation during the mixing time. The cross-peak build-up shown on the bottom left was simulated for two exchange times $1/k_1 = 50$ ms and $1/k_2 = 100$ ms. The exchange time for the curves on the bottom right is $1/k_2 = 100$ ms. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]</figcaption></figure>

## 2D Simulations: Three-Site Exchange

The results obtained for two-site exchange simulations did follow the known analytical solutions recalled in section “general analytical solutions for two-site exchange” (2, 8, 23). To date, there is no published analytical solution for the case of relaxation and exchange between more than two sites. However, this case can be modeled with Eq. [23]. There is a fundamental difference between two- and n-site exchanges. For two sites, the exchange maps are symmetric (if both dimensions are acquired with the same evolution times) as mass balance requires the forward and backward exchange to be identical. When exchange takes place between more than two sites, detailed mass balance applies to the overall exchange but not to subsets consisting of only two of the *n* sites, so that asymmetric exchange maps can arise. Figure 7 shows two simulated maps for three-site exchange, one symmetric (left) and the other one not (right). To obtain a symmetric map, all

exchange times have to have the same value, whereas they differ strongly in the nonsymmetric case.

## 2D Simulations: Noise and Baseline Artifacts

A further source of asymmetry in experimental data sets can be linked to the sensitivity of the numerical method used to approximate the inverse Laplace problem (the ILT algorithm) to noise and baseline artifacts (23). To prove this point, a constant baseline and zero-mean Gaussian white noise with a standard deviation equal to a percentage of the maximum value in the simulated 2D data sets without any noise were added to the simulated data. The corresponding exchange maps derived by 2D ILT were analyzed for the peak integrals and the peak coordinates in the $T_2-T_2$ space (Fig. 8). Figure 8(a) quantifies the effect of the noise on the peak integrals. The amplitude of the
<figure><img src="image_8.png" /><figcaption>**Figure 7** 2D maps for three-site exchange. (a) Symmetric map with $k_{12} = k_{13} = k_{23} = k_{32} = 250$ ms. (b) Asymmetric map with $k_{12} = 400$ ms, $k_{13} = 20$ ms, $k_{23} = 3000$ ms, and $k_{32} = 10$ ms. The longitudinal relaxation times were set to four times the corresponding transverse relaxation times.</figcaption></figure>

noise strongly affects the value of the integrals. As the peak integrals also depend on both the exchange rate and the longitudinal relaxation times, the accuracy with which exchange rates and longitudinal

relaxation times can be determined is affected by the noise. The peak positions can also be shifted on the $T_2-T_2'$ map in the presence of noise [Fig. 8(b)]. This effect is not negligible as, remarkably, one of the di-

<figure><img src="image_9.png" /><figcaption>**Figure 8** Effects of the noise (a and b) and baseline offset (c and d) on the peak intensities (a and c) and peak coordinates (b and d). The relative deviation of the integrals from the value in the absence of noise is given in percent. The baseline offset is specified relative to the maximum value of the original time-domain 2D data set. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9** Water saturating silica particles. (a) Transverse magnetization decay. (b) Longitudinal magnetization build-up measurements of water saturating silica particles. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]</figcaption></figure>

agonal peaks shifts away from the diagonal when the noise level is only just above 2%.

Nonzero baselines also produce similarly strong artifacts in the peak intensities (c) and the peak positions (d). As shown in Figs. 8(c,d), a baseline offset of 2% is sufficient to produce distorted exchange maps.

## IV. COMPARISON WITH EXPERIMENTS

The formalism to simulate 2D relaxation exchange maps was tested by modeling experimental maps for two-site exchange to extract relaxation and exchange parameters. Silica particles were synthesized according to the Stöber reaction (27) as this procedure is known to yield monodisperse particles, which are good candidates for simple two-site exchange when

<figure><img src="image_11.png" /><figcaption>**Figure 10** Flowchart demonstrating the procedure used to fit the experimental data. At each iteration, the simplex algorithm is used to estimate a new set of parameters minimizing the least-squares misfit of measured and experimental data points.</figcaption></figure>

saturated with water. The two relaxation sites are expected to be the bulk water in the interstitial space of the closely packed structure on the one hand, and the water on the surface and within the micropores of the silica particles on the other hand. For the particle synthesis, 60 mL of tetraethoxysilane, 405 mL of $NH_3$ at 25% in solution, and 1.5 L of ethanol were mixed together. Then, the particles were washed first with ethanol and a second time with water, and their size was determined with a scanning electronic microscope to be 500 nm with a standard deviation of 40 nm.

Figure 9 depicts the transverse (a) and longitudinal magnetization build-up (b) curves measured by the CPMG sequence and the saturation recovery method, respectively. As expected, the system exhibits two apparent transverse relaxation times $T_2$, which correspond to two different sites, but only a single apparent relaxation time $T_1$. Considering the results obtained in the 1D simulations (see section “procedure for 1D simulations”), we can conclude from the observation of a single $T_1$ that the exchange time is small compared to the spin-lattice relaxation time. In

<table><caption>**Table 4** Fit Parameters Obtained by the Simplex Algorithm</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>$M^A$</td><td>5.5</td></tr><tr><td>$M^B$</td><td>18.4</td></tr><tr><td>$T_2^A$</td><td>3.9 ms</td></tr><tr><td>$T_2^B$</td><td>103.5 ms</td></tr><tr><td>$T_1^A$</td><td>610 ms</td></tr><tr><td>$T_1^B$</td><td>630 ms</td></tr><tr><td>$k$</td><td>1/110 ms<sup>-1</sup></td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Figure 11** Experimental (a) and simulated (b) T₂–T₂ exchange spectra of silica particles with the parameters of Table 3 for a mixing time of 70 ms.</figcaption></figure>

this case, the single apparent T₁ corresponds to the nonzero eigenvalue of the **R**⁽¹⁾ + **K** matrix.

T₁–T₂ correlation and T₂–T₂ exchange experiments were performed with a homemade Halbach-magnet operating at a 21.85 MHz ¹H frequency. Using the algorithm described in section II, 2D time data set were simulated for two-site exchange with six different mixing times (0, 10, 30, 50, 70, and 90 ms). A simplex algorithm (28) was used to match the input parameters of the simulation with the experimental data for each mixing time. The data were first fitted directly in the time domain to avoid any artifact from the ILT. The fitting procedure is described in Fig. 10. The fit parameters are summarized in Table 4.

Although this is beyond the point of this study, it is interesting to note that the ratio of the observed amplitudes is 3.3. This is an order of magnitude below the expected value considering only the ratio between the interstitial space and the microporous volume. With a close packing volumetric density of 0.74 and porous volume within the particles of about 0.6% (a common

value for Stöber silica) (29), one should get a ratio of 43. There is thus a significant contribution of the surface water to the magnetization M<sub>A</sub>.

Figure 11 shows the results obtained after ILT with a mixing time of 70 ms. One (a) is obtained by inverting the experimental data set, and the other one (b) by inverting the best-fit simulated data set. The peaks of the simulated maps are broader than the experimental ones. This is interpreted as a regularization artifact of the ILT procedure. To verify the quality of the fit, the evolution of the peak integrals with the mixing time for both simulated and experimental maps is compared in Fig. 12.

Except for a discrepancy concerning the diagonal peak 1 at low mixing time, Fig. 12 exhibits satisfactory agreement between the mixing time dependence of the peak intensities observed experimentally and the ones simulated using the parameters obtained by fitting the 2D data set. Thus, fitting experimental 2D time-domain data sets with simulated ones allows one to unravel relaxation from

<figure><img src="image_13.png" /><figcaption>**Figure 12** Comparison between the integrals of diagonal (a) and cross (b) peaks obtained by simulation and experiment versus the mixing time. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 13** Experimental (a) and simulated (b) $T_1-T_2$ correlation map of water-saturated silica particles derived with the parameters of Table 3.</figcaption></figure>

exchange and obtain the true $T_2$ values, which are quite different from the apparent ones obtained by simple 1D fits [Fig. 9(a)]. Unfortunately, in the experimental system under study, it has not been possible to extract real $T_1$ values as the exchange time was too short to distinguish them.

In support of the exchange analysis, the $T_1-T_2$ correlation map was measured and compared to the simulated one obtained with the fit parameters presented in Table 4 (Fig. 13). The small peak in the experimental map at low $T_1$ (around 40 ms) is likely to be noise and does not carry any meaning as the 1D fit of Fig. 9(a) clearly indicated a single $T_1$.

## V. CONCLUSION

In this work, we have shown that even if analytical solutions are known for two sites, a numerical simulation is required to analyze exchange maps with more than two sites. The simulations allowed us to propose a source for asymmetry, which is often observed experimentally as well as a source for the shift of peaks. Moreover, we proposed a method to extract the exchange kinetic parameters by fitting the experimental 2D time-domain data sets with simulated ones. As an example, two sites with an exchange rate of $1/110 \text{ ms}^{-1}$ were identified for water in a saturated dense packing of silica particles with a narrow particle size distribution.

Although the $T_2-T_1$ exchange case has extensively been discussed, this approach can be easily extended to $T_1-T_2$ correlation maps for which one example (Fig. 13) was given. Such experiments can be understood in terms of $T_2-T_1$ experiments which start with $T_1$ encoding by saturation or inversion recovery instead of $T_2$ encoding in the evolution time, and where the mixing time vanishes.

The theoretical framework allows the study of $n$ discrete-site systems, but its use is limited in practice to the case of four or five sites. Otherwise, the number of independent parameters becomes too large. On the other hand, it provides a first step toward the study of continuous distributions of sites.

## ACKNOWLEDGMENTS

The authors express their gratitude to Prof. N. Lequeux (ESPCI) for the Stöber synthesis protocol as well as to Dr. Federico Casanova and Dr. Juan Perlo (RWTH) for truly helpful discussions. Financial support for this work was provided by Saint-Gobain Recherche (France) and TR32 (Germany).

## REFERENCES

1. Hürlimann MD, Venkataramanan L. 2002. Quantitative measurement of two-dimensional distribution functions of diffusion and relaxation in grossly inhomogeneous fields. J Magn Reson 157:31–42.

2. McDonald PJ, Mitchell J, Mulheron M, Aptaker PS, Korb JP, Monteilhet L. 2007. Two-dimensional correlation relaxometry studies of cement pastes performed using a new one-sided NMR magnet. Cem Concr Res 37:303–309.

3. Song YQ, Zielinski L, Ryu S. 2008. Two-dimensional NMR of diffusion systems. Phys Rev Lett 100:248002.

4. Washburn KE, Callaghan PT. 2006. Tracking pore to pore exchange using relaxation exchange spectroscopy. Phys Rev Lett 97:175502.

5. Callaghan PT, Arns CH, Galvosas P, Hunter MW, Qiao Y, Washburn KE. 2007. Recent Fourier and Laplace perspectives for multidimensional NMR in porous media. Magn Reson Imaging 25:441–444.
6. Callaghan PT, Godefroy S, Ryland BN. 2003. Diffusion-relaxation correlation in simple pore structures. J Magn Reson 162:320-327.

7. Chandrasekera TC, Mitchell J, Fordham EJ, Gladden LF, Johns ML. 2008. Rapid encoding of T-1 with spectral resolution in n-dimensional relaxation correlations. J Magn Reson 194:156-161.

8. McDonald PJ, Korb JP, Mitchell J, Monteilhet L. 2005. Surface relaxation and chemical exchange in hydrating cement pastes: a two-dimensional NMR relaxation study. Phys Rev E 72:011409.

9. Jeener J, Meier BH, Bachmann P, Ernst RR. 1979. Investigation of exchange processes by 2-dimensional NMR-spectroscopy. J Chem Phys 71:4546-4553.

10. Ernst RR, Bodenhausen G, Wokaun A. 1987. Principles of Nuclear Magnetic Resonance in One and Two Dimensions. Oxford: Oxford University Press.

11. Wüthrich K. 1995. NMR in Structural Biology. Singapore: World Scientific.

12. Macomber RS. 1998. A complete introduction to modern NMR spectroscopy. New York: Wiley.

13. McConnell HM. 1958. Reaction rates by nuclear magnetic resonance. J Chem Phys 28:430-431.

14. Bain AD. 2003. Chemical exchange in NMR. Prog Nucl Magn Reson Spectrosc 43:63-103.

15. Kroeker RM, Henkelman RM. 1986. Analysis of biological NMR relaxation data with continuous distributions of relaxation-times. J Magn Reson 69:218-235.

16. Lee JH, Labadie C, Springer CS, Harbison GS. 1993. Two-dimensional inverse Laplace transform NMR: altered relaxation times allow detection of exchange correlation. J Am Chem Soc 115:7761-7764.

17. Venkataramanan L, Song YQ, Hürlimann MD. 2002. Solving Fredholm integrals of the first kind with tensor product structure in 2 and 2.5 dimensions. IEEE Trans Signal Process 50:1017-1026.

18. Ernst RR, Anderson WA. 1966. Application of Fourier transform spectroscopy to magnetic resonance. Rev Sci Instrum 37:93-102.

19. Morris KF, Johnson CS. 1993. Resolution of discrete and continuous molecular-size distributions by means of diffusion-ordered 2D NMR-spectroscopy. J Am Chem Soc 115:4291-4299.

20. Callaghan PT, Manz B. 1994. Velocity exchange spectroscopy. J Magn Reson Ser A 106:260-265.

21. Sun BQ, Dunn KJ. 2005. A global inversion method for multi-dimensional NMR logging. J Magn Reson 172:152-160.

22. Hürlimann MD, Burcaw L, Song YQ. 2006. Quantitative characterization of food products by two-dimensional D-T-2 and T-1-T-2 distribution functions in a static gradient. J Colloid Interface Sci 297:303-311.

23. Monteilhet L, Korb JP, Mitchell J, McDonald PJ. 2006. Observation of exchange of micropore water in cement pastes by two-dimensional T-2-T-2 nuclear magnetic resonance relaxometry. Phys Rev E 74:061404.

24. Fleury M, Soualem J. 2009. Quantitative analysis of diffusional pore coupling from T2-store-T2 NMR experiments. J Colloid Interface Sci 336:250-259.

25. Whittall KP, Mackay AL. 1989. Quantitative interpretation of NMR relaxation data. J Magn Reson 84:134-152.

26. Straley C, Schwartz LM. 1995. Transverse Relaxation in Random Bead Packs: Comparison of Experimental Data and Numerical Simulations. Louvain, Belgium: Pergamon-Elsevier Science Ltd. pp 999-1002.

27. Green DL, Lin JS, Lam YF, Hu MZC, Schaefer DW, Harris MT. 2003. Size, volume fraction, and nucleation of Stober silica nanoparticles. J Colloid Interface Sci 266:346-358.

28. Nelder JA, Mead R. 1965. A simplex-method for function minimization. Comput J 7:308-313.

29. Szekeres M, Toth J, Dekany I. 2002. Specific surface area of Stoeber silica determined by various experimental methods. Langmuir 18:2678-2685.
BIOGRAPHIES

<figure><img src="image_15.png" /></figure>

**Maxime Van Landeghem** was born in Paris (France) in 1984. He obtained his diploma from the French graduate school ParisTech National Institute of Physics and Chemistry (ESPCI ParisTech) in 2007 and his M.Sc. in 2008. Now, he is performing his PhD under the direction of Prof. Jean-Baptiste d'Espinose de Lacaillerie (ESPCI ParisTech) and Prof. Bernhard Blümich (RWTH). His PhD is focused on the understanding of drying and hydration process in porous media. He is funded by Saint-Gobain Research.

<figure><img src="image_16.png" /></figure>

**Agnes Haber** was born in Covasna (Romania) in 1984. She obtained her B.Sc. (2006) and M.Sc. (2008) at Babes-Bolyai University in Cluj-Napoca (Romania). She is now doing her PhD at ITMC, RWTH Aachen, Germany, under the supervision of Prof. Dr. Bernhard Blümich. Her main research topic is focused on water dynamics in different porous materials by low-field NMR.

<figure><img src="image_17.png" /></figure>

**Jean-Baptiste d'Espinose de Lacaillerie** is an assistant Professor at ParisTech National Institute of Physics and Chemistry (ESPCI ParisTech) in Paris. He is a member of the Soft Matter and Engineering Laboratory (a joint laboratory between ESPCI, the University Pierre et Marie Curie and the National Centre for Scientific Research). He holds a geological engineering diploma and a M.Sc. in Geochemistry from the National Institute of Geology (ENSPG INPL) in Nancy, France (1987) and a PhD. in Chemistry from the University of Wisconsin-Milwaukee (1992). He did his postdoctoral research at the French Petroleum Institute. He has coauthored more than 70 publications in the area of catalysis, surface chemistry of oxides in aqueous solutions, and solid-state NMR. His current research interests focus on the physics and chemistry of organomineral interactions at the nanoscale, the chemistry of cements substitutes for low carbon emission, and stray-field imaging of drying porous media.

<figure><img src="image_18.png" /></figure>

**Bernhard Blümich** is Chair of Macromolecular Chemistry at RWTH Aachen University. His research interests are the development and applications of NMR for analysis of materials and processes. Starting with NMR as an undergraduate at the Technical University Berlin, he has since contributed to the methodological development of NMR in the areas of multidimensional stochastic NMR spectroscopy, solid-state NMR of polymers, NMR imaging of materials, and mobile NMR. He has published over 300 articles, edited six books, and written two monographs. In 2006, he was Miller Professor at the University of California at Berkeley, received the Ampere Price 2007 for the development of the NMR-MOUSE and its applications to various fields of science and technology, and was nominated Fellow of the International Society of Magnetic Resonance in 2008.