Dipolar Relaxation of Water Protons in the Vicinity of a Collagen-like Peptide

Quantitative magnetic resonance imaging is one of the few available methods for noninvasive diagnosis of degenerative changes in articular cartilage. The clinical use of the imaging data is limited by the lack of a clear association between structural changes at the molecular level and the measured magnetic relaxation times. In anisotropic, collagen-containing tissues, such as articular cartilage, the orientation dependency of nuclear magnetic relaxation can obscure the content of the images. Conversely, if the molecular origin of the phenomenon would be better understood, it would provide opportunities for diagnostics as well as treatment planning of degenerative changes in these tissues. We study the magnitude and orientation dependence of the nuclear magnetic relaxation due to dipole–dipole coupling of water protons in anisotropic, collagenous structures. The water–collagen interactions are modeled with molecular dynamics simulations of a small collagen-like peptide dissolved in water. We find that in the vicinity of the collagen-like peptide, the dipolar relaxation of water hydrogen nuclei is anisotropic, which can result in orientation-dependent relaxation times if the water remains close to the peptide. However, the orientation-dependency of the relaxation is different from the commonly observed magic-angle phenomenon in articular cartilage MRI.


INTRODUCTION
Articular cartilage (AC) is the soft tissue that covers the bones at joint surfaces. It reduces friction and provides shock absorption in the joint motion. Magnetic resonance imaging (MRI) is a noninvasive method well suited to study soft tissues and holds a promise for being an efficient diagnostic tool for early degenerative changes in articular cartilage. However, the association of the molecular-level structure to quantitative MRI parameters, such as the transverse relaxation time T 2 , remains an open problem.
The most abundant macromolecular component in articular cartilage is type II collagen, which amounts to 20−40% of the wet weight of cartilage. Collagen consists of tropocollagens that form prototypic fibrils with a thickness of 18 ± 5 nm, which then bundle up to larger fibers with diameters typically below 100 nm in knee articular cartilage. 1 The tropocollagen unit consists of three polypeptide strands, which form a triplehelical structure. In addition to collagen, the other important macromolecular component is proteoglycans, which are embedded in and interact with the collagen matrix. The majority of the wet weight of articular cartilage is water, which surrounds and penetrates the fibers/fibrils. 2 1 H MRI of articular cartilage reflects the interactions of water with the macromolecular structures it surrounds. Articular cartilage can be divided into three zones by the orientation of collagen fibrils. 2 In the superficial zone, fibers are parallel to the surface. Deeper, in the radial zone, the fibrils are on average aligned perpendicular to the bone-cartilage interface. Between the superficial and radial zones is the transitional zone, in which the fibril orientation is more or less random. Changes in collagen orientation, especially in the superficial zone, have been linked to early osteoarthritis. 3 The anisotropy of relaxation times in tissues containing aligned networks of collagen has been studied experimentally in cartilage 4−7 and tendon, 4,8 with results mostly pointing to either weak or negligible anisotropy in the longitudinal relaxation time T 1 . For T 2 -weighted images and T 2 maps, the story is different: If one rotates a small cartilage sample around an axis perpendicular to the bone-cartilage interface, the T 2weighted image of the radial zone is largely unaffected. Rotations around the other axes produce systematic variations in T 2 -weighted MRI images, and the corresponding relaxation time maps. 6,7,9 This phenomenon is known as the magic angle effect 10 and in collagenous tissues it is usually attributed to nonaveraged dipolar couplings of water bound to the collagen. 8,11,12 The orientational symmetry of T 2 in the radial zone suggests that the average macromolecular environment experienced by the hydrogen atoms of water can be assumed to be cylindrically symmetric.
Competing explanations for the relaxation anisotropy include the formation of strong water bridges 11 and anisotropic, water-filled cavities inside the tissue. 13,14 In both models, the anisotropy in relaxation is due to residual 1 H− 1 1H dipole−dipole couplings, although they arise for different reasons. The role of proton exchange in relaxation is brought forward, especially in studies of the longitudinal relaxation time in the presence of a continuous-wave spin-lock (T 1ρ ). 15 However, the magic-angle effect is seen also in T 1ρ of articular cartilage, 9 and a comprehensive theoretical picture of how proton exchange would create the orientation-dependence is lacking.
Most of the earlier computational work on collagen has focused on the structure and biochemistry of collagen itself. 16 Anisotropic rotational motion of water near collagen has been analyzed previously by numerical and analytical models, 12 as well as molecular dynamics (MD) simulations focusing on the so-called water bridges. 17,18 Water diffusion in a collagen network has been computationally modeled with coarsegrained Monte Carlo simulations of Brownian dynamics 19 and Langevin dynamics. 20 Although the effect of proteins and other solutes on water dynamics is an extensively studied subject, 21−23 these studies do not consider macromolecules in a fixed or restricted orientation condition as they appear in anisotropic tissues, such as articular cartilage. Furthermore, little, if any attention is given to the effects which could be seen in 1 H NMR/MRI due to the anisotropic motion of water in these systems.
In this paper we study the sources of 1 H relaxation anisotropy with MD simulations of a small collagen-like peptide dissolved in water. With our molecular model we aim to estimate how the dipolar relaxation would appear in a uniaxially anisotropic network of collagen molecules, such as the radial zone of articular cartilage. In the present work we try to avoid assumptions about the motional characteristics or specific binding sites of the water molecules and compute relaxation rates directly from the simulated data. Our approach represents a first step toward more comprehensive and realistic models of water in a collagenous environment. The validity of the methodology for computing the dipolar relaxation from MD simulations has previously been demonstrated for pure water. 24 In the Supporting Information (SI) we derive the necessary extensions to the well-known Redfield relaxation theory 25−27 of dipolar relaxation in isotropic liquids to handle the residual dipole−dipole couplings (RDCs) and to study the orientation dependence of relaxation with respect to the direction of the main magnetic field. Unlike in the seminal work by Woessner, 28,29 we start our derivation from first-principles, without the assumption of how the relaxation rates depend on the spectral densities of molecular motion. We also do not assume any overall motion that would average out the RDC. Similar developments, albeit with a slightly different formalism have been given, for example, in ref 30. With the combination of theoretical work and simulations, we clarify the role of water rotation and dipolar interaction in the 1 H relaxation times in quantitative MRI of articular cartilage and other tissues containing aligned networks of collagen.

THEORY
Here we estimate the evolution of the longitudinal and transverse magnetizations in the presence of time-dependent and residual dipole couplings. The equations are derived for a pair of identical spin-1/2 nuclei. We present here mainly the results. The complete derivation can be found in the SI. The spin Hamiltonian can be split to time-independent and timedependent parts H 0 and H 1 (t) as where H Z = ω 0 (I z + S z ) and H DD refer to Zeeman and dipole− dipole interactions, respectively, for the two spins I and S. The time-average ⟨H DD (t)⟩ t is the RDC. The dipole−dipole coupling can be decomposed into spatial functions F m and spin-operators V m so that where r t = r(t), is the length of the internuclear vector and θ t = θ(t) and ϕ t = ϕ(t) describe the orientation of the vector with respect to the external magnetic field. In SI units C = −2.93083474 × 10 −24 m 3 /s. The spin operators V m with |m| = 0, 1, and 2 correspond to in-phase magnetization with longitudinal two-spin order, antiphase magnetizations, and double-quantum coherences, respectively: To aid in interpretation of the molecular motion, the spherical harmonics can be written in Cartesian coordinates, here for brevity only in the hemisphere 0 < θ < π/2: We see that if r is constant as in intramolecular coupling with a rigid water molecule model, the correlation time with m = 0 is dependent only on the z coordinate, m = ±2 is dependent only on motion in the x, y-plane, and m = ±1 represents a mixed motion involving movement in-plane and also along the zdirection.
We assume that the motion of internuclear vectors r has on average uniaxial symmetry around the long axis of the collagenlike peptide. We analyze the validity of this assumption in section S7 of the SI. A coordinate system with one axis, denoted z′-axis, along the symmetry axis is taken as the principal axis system (PAS) of the RDC. Treating the dipolar coupling as a small time-dependent perturbation with Redfield theory (see sections S1 and S2 in the SI) gives the evolution of the longitudinal magnetization when the correlation time(s) τ c of the dipolar coupling are much shorter than the inverse of the secular RDC. The spectral densities in any coordinate system tilted with respect to the PAS by an angle β can be expressed with the spectral densities in the PAS, J′ nn (ω), using the rotation formula (section S5 in the SI) Here, β is the angle between the symmetry axis and the main magnetic field. d m,n (2) (β) are elements of the Wigner small-d matrix. From here on, quantities denoted with a prime (′) are computed in the coordinates in which the system is uniaxially symmetric around the z′-axis. The spectral densities in the PAS are The time-correlation functions (TCFs) G nn ′ (|τ|) are computed in the PAS of the RDC, The overline denotes ensemble average. θ′ is the angle between the axis of symmetry and the internuclear vector. ϕ′ is the other spherical coordinate in the plane perpendicular to the symmetry axis. Again, we use the assumption of a cylindrically symmetric macromolecular environment for the water molecules and choose the x′ and y′ axes in this plane arbitrarily. The subscripts t and t + τ denote the respective time points. For transverse magnetization we find (see sections S3 and S4 in the SI) the time-evolution of the rotating frame operator I x r + S x r to be of the form In addition to the spectral densities encountered earlier, eq 14 contains the effect of the secular RDC in the term Δ = 3d IS , where d IS the secular dipole−dipole coupling as defined in ref 31.
The transverse magnetization is a product of an exponential decay and two oscillating terms. In section 5.5 we see that in our case the oscillation by cos {−1/8 [L 1,1 (ω 0 ) + 2L 2,2 (2ω 0 )] t} is much slower than the one caused by the term containing the secular RDC. It is also much slower than the relaxation and can be neglected to a good approximation.
If the RDCs are nonvanishing, the form of the evolution of the transverse magnetization depends on the sign of 4Δ 2 − 0,0 1,1 0 2,2 0 If also the RDC vanishes (Δ = 0), the relaxation rate of the transverse magnetization becomes the familiar 25,26 3. MODELING 3.1. The Molecular Model. Here we chose the trimer of the polypeptide (Pro-Hyp-Gly) 4 -Glu-Lys-Gly-(Pro-Hyp-Gly) 5 with PDB ID: 1QSU 32 to represent a partial sequence of the triple-helical structure of a tropocollagen of type II collagen. The length of the peptide in its initial state is 8.8 nm. The typical length and diameter for a type-II collagen monomer are 1.5 and 300 nm, respectively. 1 In articular cartilage, the diameter of a collagen fiber can be up to 200 nm. For atomistic simulation these structures are huge, which makes their simulations challenging. We can, however, quite safely assume that the most important collagen−water interactions affecting water dynamics are short-ranged. In effect, the individual water molecules mostly experience the local structure of the collagen. This should allow us to model the collagen with small representative samples of the fibrillous structure. The same The Journal of Physical Chemistry B pubs.acs.org/JPCB Article may not apply for the intermolecular contribution to the residual dipolar couplings. We will discuss this limitation further in section 4.1. A weak long-ranged (∝ r −3 ) aligning force felt by the water molecules can be produced by the electric field of the protein. 21 This could be reflected in the intramolecular 1 H− 1 H dipole−dipole couplings in water. The long-range contribution from electrostatic forces is accounted by the particle mesh Ewald method as explained in section 3.2.
Since the long-range forces are computed with a periodic approximation in the PME method, the size and shape of the simulation box and the proportion of water will affect the results. We did not try to separate the contributions of shortand long-range forces to the water alignment. In simulations, the peptide was solvated in TIP4P/2005 water, which is a rigid, nonpolarizable, planar, four-site model. 33 We assume that bond vibrations are mostly uncoupled from rotational motion and that the vibrations themselves are too fast to cause relaxation. Among the commonly used rigid, nonpolarizable water models, TIP4P/ 2005 has been shown to best reproduce the experimental rotational correlation times in simulations. 34 Relatedly, it is also well suited for simulating 1 H nuclear spin relaxation. 24 The initial coordinates as well as the force field parameters were created with xLeap, which is part of the AmberTools15 package. 35 The Amber force field ff14SB 35 with rigid bonds was used to model the atomic interactions. A cutoff radius of 12 Å was used for nonbonded forces, with a switching function applied at 10 Å. Before the solvation, xLeap was used to add hydrogen atoms and missing heavy atoms to the peptide. This left the C-termini glycine residues deprotonated and the Ntermini proline residues protonated.
3.2. Simulations. The MD simulations were performed with NAMD. 36 Periodic boundaries were applied in all directions. For the bulk water simulation, the cubical simulation box contained 1926 water molecules and the initial box side length was 43 Å. The system with the triple-helical peptide contained 18 393 water molecules. The simulation box side lengths were initially L x , L y , L z = 86.7 Å, 60.5 Å, 128.6 Å, respectively.
The simulations of both systems consisted of 1000 steps of potential energy minimization, followed by ≈18 ns of constanttemperature, constant-pressure (NPT) ensemble simulation at T = 300 K and P = 1 atm. The Langevin thermostat with 1/ps damping factor was used to control temperature. The Langevin-piston-inspired barostat present in NAMD 36 was used with a 100 fs oscillation period and 50 fs decay. A time step of 1 fs was used for the integrator, and nonbonded forces were computed at each step. Long-range electrostatic forces were updated every second step with the particle mesh Ewald method. Neighbor lists with 14 Å lookup radius were used for nonbonded interactions. The lists were updated every tenth time step. The trajectory of the molecules was recorded every 100 fs. First 3 ns was discarded as equilibration, and the rest of the trajectory was used for analysis. The input files to reproduce the simulations are available for download at the Zenodo archive (10.5281/zenodo.6330600). The trajectories and log files of the simulations are available from the corresponding authors by request.

ANALYSIS OF THE TIME-CORRELATION FUNCTIONS OF DIPOLE−DIPOLE COUPLINGS
We assume that on average the dipolar couplings experienced by the water molecules are cylindrically symmetric and the symmetry axis is defined by the eigenvector corresponding to the smallest eigenvalue of the inertia tensor of the peptide. The correlation functions in eq 12 were computed in the average PAS of the inertia tensor of the peptide. Then the functions were fitted with a triexponential, The multiexponential correlation function leads to a spectral density This form can easily be used to compute the relaxation rates, once the parameters A mk and τ mk have been obtained from the fit.
Since the trajectory sampling is every 100 fs, the details of the femtosecond scale libration of the water molecules do not show up in the correlation functions. We consider the contribution of this motion to relaxation as small. If one monoexponential function in eq 20 would be dedicated to describing the libration, its amplitude A mk would be small and its correlation time τ mk very short. Other exponential terms corresponding to, for example, the overall rotation of the molecule, have both larger amplitude and longer correlation time and they dominate the spectral density of eq 20. The libration can and will affect the time-evolution of the magnetization by decreasing the RDCs.
Three exponentials were found sufficient to capture well the characteristics of the correlation functions in most cases studied here, with one exception (see section 5.7). To estimate the error made in fitting the TCFs we compared the integrals of the TCFs and their respective fits, which commonly gave a difference of 2% or less. The number of exponentials in the fit is not relevant as long as it does not substantially affect the accuracy of the fit. We do not claim that the individual τ mk represent for example different characteristic motions. Instead, we compute the integrated correlation times to analyze the molecular motion in the PAS of the inertia tensor. These integrated correlation times are computed in section 5.2 and their connection to relaxation anisotropy is discussed in section 6 via the Cartesian representation of spherical harmonics [eq 8]. It would also be possible to proceed from the correlation functions to the relaxation rates by numerically computing the Fourier transformations of the correlation functions as has been done in the case of 129 Xe relaxation due to chemical shift anisotropy. 37 Our strategy allows the analysis of the integrated correlation times and to draw the connection from them to the relaxation rates. The correlation functions and RDCs appearing in the relaxation eqs 9, 14, and 18 were computed in the principal The Journal of Physical Chemistry B pubs.acs.org/JPCB Article axis system of the inertia of the peptide. The assumed cylindrical symmetry implies that there is only one angle needed for the transformation of results for different peptide orientations with respect to the main magnetic field. This angle β is visualized in Figure 1.

Water Selections.
The analysis of relaxation in the vicinity of the peptide is separated into contributions from water in different geometrical regions around the peptide, as well as contributions from different types of coupling. Intraand intermolecular couplings in water are referred to as intra and inter, respectively. The couplings between hydrogen nuclei in water and in the peptide are denoted with pept. H. The sum of intra, inter, and pept. is referred to as all. The different contributions are additive, provided that the cross-correlations between the corresponding dipolar couplings are negligible. Water molecules which have any atom of the collagen-like peptide closer than 3.5 Å are denoted with 1st h.l. (h.l. = hydration layer) (Figure 2). Similarly, 2nd h.l. contains the water molecules which have any atom of the peptide at a distance r, r < 7.0 Å but not at r < 3.5 Å. Finally, water contains all water molecules in the simulation with the peptide, for example, water pept. H refers to 1 H− 1 H couplings between all water molecules and the peptide. The criterion for the hydration layer, r < 3.5 Å, contains the first two peaks of the water−peptide O−H and H−O radial distribution functions (RDFs), as well as the first peaks of the O−N and O−O RDFs ( Figure S1, in the SI). Our definition for the second hydration layer is chosen on RDFs to fully contain the second hydration layer at the cost of also containing contributions of higher hydration layers. The definition of the hydration layers is rather coarse (see section S6 in the SI) but adequate for the purpose of seeing the effect of the vicinity of the peptide to the molecular motion of water. The results labeled with bulk are computed from the bulk water simulation, including all water molecules therein.
Spherical cutoff radii of 10 and 20 Å were used for the intermolecular couplings in the bulk and peptide simulations, respectively. In the system with the peptide, the simulation cell is anisotropic. Computing the intermolecular couplings between water molecules without cutoff would bias the RDCs toward the directions in which the simulation cell is the longest. Compared to the RDCs, the TCFs decay quickly with distance and the long-range effects are much less of concern. Although dipolar couplings decay with distance as r −3 , their sum does not. So, as opposed to the correlation functions, the RDCs are long-ranged and we can not estimate them reliably from the simulation except in the case of intramolecular couplings. This limits the use of eq 14 for our analysis.
By separately studying the relaxation rates of protons in water molecules at different locations relative to the peptide, we aim to see how the vicinity of the macromolecule affects the water proton relaxation in a fibrillar environment such as articular cartilage. Average distance between collagen type II monomers in a fibril is 13 Å, 38 which means the water inside fibrils is confined very close to the macromolecules until it eventually exchanges with the bulk water outside. The space accessible for water molecules inside the fibrils corresponds roughly to the combined first and second hydration layers in our study.
4.2. Analysis Based on Initial Location of Water Molecules. Two types of analysis was employed for the computation of the TCFs. In the first one, which we later refer to as "unrestricted" analysis, the trajectory given by the simulations was sliced into 100 ps blocks. The correlation functions in eq 12 were computed in the average PAS of the inertia tensor of the peptide in each of these blocks. The water selections used in this type of analysis are static and defined based on the location of the water molecule at the beginning of the 100 ps trajectory segment. This simplified definition of hydration layers allows for a more thorough statistical analysis of the results (i.e., bootstrapping) compared to the more stringent definition used in the restricted analysis below. Also, together with the restricted analysis, it can give an assessment  The Journal of Physical Chemistry B pubs.acs.org/JPCB Article of the impact of residence times in the hydration layers on TCFs. In practice the correlation function in this case becomes where t 0 is the beginning of the block and the latter term represents the time-independent baseline. The water molecules in the hydration layers of the peptide eventually escape the space corresponding to their assignment. During the first 100 ps, more than 90% of the water molecules initially in the first hydration layer have left it at least once and after the first nanosecond all water molecules have left the first hydration layer. Since we compute the correlation functions as averages over 100 ps blocks, a water molecule initially in the first hydration layer has most likely spent time in the second hydration layer and possibly beyond. For these reasons, the ⟩| of the correlation functions in this type of analysis are always computed for all water molecules. For the same reasons the secular RDCs in Δ were computed for all water molecules. The starting point for the TCF calculation, t 0 , is not iterated through the 100 ps block. Therefore, the TCF is essentially independent of block length. In section 5.2 it is seen that the correlation times of the dipole−dipole couplings computed using the trajectory segments are mostly much shorter than the length of the segment and the error due to segmenting should be relatively small.
Where shown, the confidence in the computed correlation times and relaxation rates, etc. was evaluated by bootstrapping. 39 The 100 ps trajectory segments used in the analysis were taken as independent, giving 150 segments in the 1QSU system. From each segment, correlation functions and RDCs were computed. This combination of quantities computed from one segment was treated as a sample. Out of all samples, 150 were randomly selected with replacement, and the correlation functions and RDCs of these samples were averaged to create a bootstrap average. This resampling scheme was repeated 1000 times to give an equal number of bootstrap averages. The desired quantities, such as spectral densities or relaxation rates at specific β were computed for each of the 1000 averages. The 2.5th and 97.5th percentile of the distribution of the bootstrap averages were used as the limits of the 95% CI. The mean of the distribution is reported as the mean value of the desired quantity. We acknowledge that the 100 ps trajectory segments are not completely independent, that is, they are correlated due to the much longer time scale of the peptide motions as compared to water molecule motions. Bootstrapping was not used to compute the PAS of the inertia tensor.

Analysis Restricted to the First Hydration Layer.
To indirectly study the effect of water molecules escaping the first hydration layer we performed also a more limited analysis in which the trajectory was used as a whole and the correlation functions were computed only during the times the water molecule stays in the first hydration layer. Results from this "restricted" analysis have the additional label restr. Here, the correlation functions were computed as where both t and t + τ iterate through the entire trajectory, but only those intervals [t, t + τ] during which the molecule stays in the first hydration layer were included in the analysis. The procedure has the effect of giving a higher weight to those water molecules which stay longer times in the hydration layer. We can consider the analyses restricted to the first hydration layer and the unrestricted analysis described in section 4.2 as two limiting cases in what can be observed with our model. Error estimates are not given for the quantities computed from the restricted analysis, since the bootstrapping methodology used for the trajectory segments could not be applied in the case when the trajectory was used as a whole.

RESULTS
5.1. Residual Dipolar Couplings. If one looks at the intramolecular RDC only and restricts the analysis to the first hydration layer so that the hydration layer is computed for each snapshot separately, we get an estimate for the RDC of water molecules that would stay close to the peptide. In this case the secular dipolar coupling is d IS = 2.66 kHz. Similar computation for the second hydration layer yields 347 Hz and for all water molecules 85 Hz. From these figures, it is evident that the RDC decays fast when moving away from the peptide, but due to the fairly strong coupling close to the peptide, the averaged interaction is non-negligible even when all water molecules are taken into account.
5.2. Correlation Times of Dipole−Dipole Couplings. Integrated correlation times [eq 21] of dipole−dipole couplings were computed from the simulations with the two different analysis methods: the unrestricted (section 4.2), in which the molecules can escape the volume of their initial assignment (e.g., first hydration layer) and the restricted (section 4.3), in which the molecules were considered in the analysis only during the time they stay in the volume of their initial assignment. Table 1 shows that the integrated correlation times [eq 21] for intra-and intermolecular 1 H− 1 H couplings between water hydrogen nuclei in the unrestricted analysis are always below 10 ps. The 1 H− 1 H couplings between water and the peptide have longer correlation times, up to 36 ps. The relatively long correlation times of the water−peptide couplings also raise the τ m ′ of the combined couplings in the first hydration layer (1st h.l. all) The Journal of Physical Chemistry B pubs.acs.org/JPCB Article close to 10 ps. The τ m ′ are mostly independent of m, which implies isotropic relative motion of the nuclei, except in the first hydration layer, where τ 2 ′ = 9.7 −0.5 +0.9 ps is slightly longer than τ 0 ′ = 8.1 ± 0.4 ps and τ 1 ′ = 7.9 ± 0.3 ps for the intramolecular couplings in water. For the intermolecular couplings between hydrogen nuclei in water and the peptide τ 2 ′ = 27.0 −0.5 +0.6 ps is shorter than τ 0 ′ = 31.3 −0.7 +0.9 ps and τ 1 ′ = 32.5 −0.8 +0.9 ps. For the couplings between hydrogen nuclei in all water molecules with the hydrogen nuclei in the peptide τ 1 ′ = 34.8 −0.8 +1.1 ps is longer than τ 0 ′ = 30.0 −0.5 +0.6 ps and τ 2 ′ = 29.9 ± 0.6 ps. The combined effect of the different couplings in the first hydration layer is that τ 0 ′ and τ 1 ′ are the same (9.4 ± 0.3 ps) and τ 2 ′ = 10.1 −0.3 +0.4 is slightly longer.
When the analysis of the couplings is restricted strictly to the water molecules in the first hydration layer of the peptide and thus the effect of water molecules escaping the vicinity of the peptide is eliminated, the correlation times for intramolecular couplings are much longer and range from τ 1 ′ = 39.3 ps to τ 2 ′ = 150 ps (Table 1). Remarkably, τ 2 ′ is much longer than τ 0 ′ and τ 1 ′. The intermolecular couplings between water are less affected and the τ m ′ with different m are very close to each other at around 17 ps. The correlation times for the water−peptide couplings are the longest ranging from τ 2 ′ = 377 ps to τ 0 ′ = 468 ps, roughly an order of magnitude longer than when only the initial position of the water molecules was considered in the analysis. Similarly to the nonrestricted analysis, τ 2 ′ is shorter than τ 0 ′ or τ 1 ′ for the water−peptide couplings. MHz were computed from the fitted correlation functions as explained in section 4. The Larmor frequency used corresponds to field strength B 0 = 9.4 T. The spectral densities show how the different molecular motions transfer to the orientation-dependency of the relaxation rates. From Figure 3 we see that for most of the water selections, the spectral densities are relatively independent of the peptide orientation. The magnitudes of J m,m with different m are also close to each other for the same contribution. These findings are confirmed in Tables S1 and S2 of section S8 in the SI, where we take a more detailed look at the spectral densities. The notable exception to the mostly isotropic behavior are the J m,m of the intramolecular couplings in water in the first hydration layer of the peptide (Figure 3). Here, each J m,m with different m has a different orientation-dependency. The J 0,0 (0) and J 1,1 (ω 0 ) have maxima at β = 90°, whereas J 2,2 (2ω 0 ) has a minimum at the same orientation. Since all spectral densities contribute to the relaxation rates with a positive sign, the anisotropies partially cancel each other.

Longitudinal Relaxation.
Longitudinal relaxation rates R 1 at ν 0 = 400 MHz were computed from the simulated trajectory using equation 9 and the methods presented in section 4. The results are listed in Table 2, from which one can see that the fastest individual relaxation rate, R 1 = 0.63 ± 0.02 s −1 is found for the intramolecular couplings in water in the first hydration layer of the peptide. The total relaxation rate in the first hydration layer due to inter-and intramolecular dipolar couplings in water, as well as intermolecular couplings between water and the peptide is R 1 = 1.12 ± 0.02 s −1 . The intramolecular couplings in water account for 56% of the R 1 in the first hydration layer. The intermolecular couplings in water, as well as between water and the peptide, contribute the remaining 26% and 18%, respectively. In the second hydration layer, the R 1 values due to inter-and intramolecular couplings of water are almost equal to each other and the contribution from the water−peptide couplings is only 6% from the total R 1 of 0.463 s −1 . When all water molecules are taken into account the total relaxation rate is R 1 = 0.322 s −1 . The intra-and intermolecular couplings in water contribute 57.3% and 40.5% to this and the water−peptide couplings provide the remaining 2.2%.
From Figure 4 we see that most of the relaxation rates are isotropic. The notable exception is the R 1 for intramolecular dipole−dipole couplings for the water molecules in the first hydration layer of the peptide. There, a clear minimum is observed when the peptide is oriented perpendicular to B 0 . From the spectral densities in Figure 3 we can see that this minimum is mostly due to J 2,2 (2ω 0 ). The anisotropic part of R 1 is 10 ± 2% of the whole relaxation rate (Table 2). Already in the second hydration layer this anisotropy vanishes (drops below 2%) and becomes insignificant when all water molecules are taken in to account. 5.5. Evolution of the Transverse Magnetization. Table  2 shows that R 2 behaves mostly similar to R 1 and therefore we will not repeat all the findings made in section 5.4. The values of the longitudinal and transverse relaxation rates differ mainly in the intramolecular couplings of water in the first hydration layer, where R 2 = 0.57 ± 0.01 s −1 is slightly slower than R 1 = This difference is translated to the overall relaxation rates in the first hydration layer. Small differences between the two relaxation rates appear also in the water− peptide couplings, but their overall effect in the relaxation is small. Notably, R 2 for the combined couplings in the first hydration layer is less anisotropic than R 1 (Table 2), and any anisotropy in the former is masked by the errors resulting from fitting the TCFs. RDCs were computed for all water molecules, similarly as in the case of baselines of the TCFs (section 4.1), and components of eq 14 were evaluated individually ( Table 3). The decay rate, R ⊥ [eq 17] in Figure 5 is almost invariable with orientation, except for the intramolecular couplings in the first hydration layer, where the rate has a weak maximum at 90°. Overall, R ⊥ is isotropic in the first hydration layer.
The term J 4 ( ) 1 4 2 1,1 0 2 ω Δ − defining the oscillation frequency in eq 14 is dominated by the secular RDC (Δ/3) and therefore its orientation dependency follows quite closely a |(3 cos 2 β − 1)| form as seen in Figure 5. The term inside the square root is positive, except perhaps very close to the magic angles 54.74°and 125.26°, where RDC vanishes. Due to partial cancellation of different contributions, the RDC from all couplings combined is smaller in magnitude than the individual contributions from water−peptide or intermolecular couplings in water. Long-range contributions to the intermolecular couplings were not evaluated, and therefore any values depending on them, such as the combined couplings in Figure 5 should be taken merely as suggestive. The other oscillating term in eq 14 has angular frequency −0.021 s −1 < 1 / 8 [L 1,1 (ω 0 ) + 2L 2,2 (2 ω 0 )] < 0.021 s −1 , which is very low. This means that we can treat the corresponding cosine term in eq 14 as equal to 1, as anticipated in section 2.
5.6. Orientation-Dependency of Relaxation Times. Relaxation times T 1,2 can be calculated from the corresponding relaxation rates R 1,2 . For bulk water, T 1 ≈ T 2 = 3.4 s. For the water with the peptide we saw in sections 5.4 and 5.5 that only the total R 1 in the first hydration layer was slightly anisotropic when the uncertainties in determining the relaxation rates were taken into account. The corresponding maximum and minimum values are T 1 (90°) = 0.93 ± 0.01 s and T 1 (0°) = 0.89 −0.02 +0.01 s, respectively. The transverse relaxation time in the first hydration layer is T 2 = 0.93 ± 0.01 s. In the second hydration layer the relaxation times agree within errors: T 1 = 2.16 ± 0.03 s and T 2 = 2.19 ± 0.02 s. For all water molecules in the peptide simulation, longitudinal and transverse relaxation times are equal within errors. T 1 = 3.10 s, which is 10% faster than the simulated longitudinal transverse relaxation times for bulk water, T 1 = 3.42 s. 5.7. Relaxation in the Analysis Restricted to the First Hydration Layer. Relaxation rates computed with the analysis restricted to the first hydration layer (section 4.3) are generally faster than in the unrestricted analysis ( Table 2). The fastest relaxation rates are found for the intramolecular couplings. Here R 1 and R 2 share the same maximum value 2.40 s −1 , which occurs for R 1 at β = 90°and for R 2 at β = 0°( Figure  6). The minimum values for the two relaxation rates are different with R 1 minimum 1.73 s −1 appearing at β = 0°and R 2 minimum 2.07 s −1 seen at β = 90°. The relaxation rates due to water−peptide 1 H− 1 H couplings result in relaxation rates of slightly smaller magnitude, but stronger anisotropy ( Table 2). The maximum and minimum R 1 due to the water−peptide couplings are 1.95 s −1 and 0.637 s −1 , at β = 90°and 0°, respectively ( Figure 6). The R 2 due to the water−peptide couplings is less anisotropic and the anisotropy is inverted as compared to R 1 . Maximum and minimum R 2 are 1.95 s −1 and 1.29 s −1 , at β = 0°and 90°, respectively. The relaxation rates due to intermolecular 1 H− 1 H couplings between water molecules are, in essence, isotropic with the differences between the orientations falling below 1%. The longitudinal and transverse relaxation rates are also practically equal to each other, with R 1 = R 2 = 0.617 s −1 at their respective maxima.
The relaxation times computed by summing up the relaxation rates from the restricted analysis give both T 1 and T 2 a minimum of 0.201 s (Figure 7). The shortest T 1 appears at β = 90°, whereas the shortest T 2 is seen at β = 0°. The Table 2. Relaxation Rates and the Fractions of Their Anisotropic Parts g(R i ) = 1 − R i,min /R i,max Computed at 400 For the peptide system, R 1 and R 2 were computed in the PAS of the inertia tensor (β = 0). For the bulk system, axes of the coordinate system were aligned along the vertices of the simulation box. The g(R i ) were computed from the minimum and maximum values of R i with respect to the angle β between the long axis of the peptide and B 0 . The errors represent 95% confidence intervals computed with the bootstrap method. See section 4.1 for an explanation of the labels.
The Journal of Physical Chemistry B pubs.acs.org/JPCB Article maxima of the two relaxation rates appear at orientations with 90-degree differences to the minima. At their respective maxima, T 1 = 0.336 s is 1.3 times longer than T 2 = 0.252 s.  Errors represent 95% confidence intervals computed with the bootstrap method. See section 4.1 for an explanation of the labels. The oscillation factor 1 / 4 (4Δ 2 − J 1,1 (ω 0 ) 2 ) 1/2 is reported only for all water molecules, since the same secular RDC values (Δ/3) were used for the analysis of the hydration layers. See section 4.1 for discussion.

DISCUSSION
We have computed relaxation of water hydrogen nuclei due to 1 H− 1 H dipolar coupling from a MD simulation of a collagenlike peptide dissolved in water. Results from a bulk water simulation were used as a reference. Redfield relaxation theory including RDC was derived for a pair of spin-1 / 2 nuclei (section 2 and SI sections S1−S5). In the system with the peptide, our two analysis strategies, the unrestricted and the restricted, could be considered as estimates of the opposite extremes. In the former, the water is fairly free to self-diffuse, and in the latter case it is confined to the macromolecular structure.
Only weak anisotropy was found in the relaxation times computed with the unrestricted analysis. The correlation times of the dipole−dipole couplings between water hydrogens were below 10 ps even in the first hydration layer. The correlation times of the couplings between hydrogens in water and the peptide were slightly longer, but still below 40 ps. The relaxation times T 1 and T 2 were 0.9 s in the first hydration layer and approached bulk water relaxation times further away from the peptide. In the analysis restricted to the first hydration layer, the simulated correlation times of the intramolecular dipole−dipole couplings varied from τ 1 ′ = 39.3 ps to τ 2 ′ = 150 ps (Table 1), reflecting significant anisotropy. For the water−peptide 1 H− 1 H couplings, the longest correlation time was on the order of 470 ps, roughly 2 orders of magnitude longer than in bulk water. The resulting relaxation times varied from 201 to 336 ms, depending on peptide orientation.
The anisotropy in the correlation times can be analyzed by looking at the spherical harmonics in Cartesian coordinates [eq 8]. The z-axis corresponds to the vector along the long axis of the peptide (Figure 1) with the x, y-plane perpendicular to this vector. For the intramolecular couplings in the first hydration layer, τ 1 ′ as the shortest correlation time would fit together with a diffusion in a cone model, in which one of the hydrogen atoms in the water molecule is bonded to the peptide via a hydrogen bond and the other hydrogen is moving more freely. This would make the mixed motion, involving movement along the z′ axis as well as in the x′, y′ plane, the fastest out of the three options. The anisotropy in the intermolecular couplings with the peptide is harder to interpret, as translational motion is also included. However, one could relate the shortest correlation time τ 2 ′ to the fact that for a water molecule adjacent to the peptide, the x, y-plane provides a way out of the peptide's proximity to diffuse into the bulk water.
6.1. Comparison of the Unrestricted and Restricted Analysis. The relaxation rates computed from the restricted analysis are generally faster than in the unrestricted analysis. The differences in the relaxation rates between the two analysis methods are the smallest in the intermolecular couplings between water hydrogens. Even there, the restricted analysis results in R 1 and R 2 twice as fast as in the unrestricted analysis. For intramolecular couplings, the relaxation rates are three to four times faster than in the unrestricted analysis, in which only the initial positions of the water molecules were taken into account. The most dramatic difference is seen in the relaxation rates resulting from the water−peptide 1 H− 1 H couplings, where R 2 is almost an order of magnitude faster in the restricted analysis as compared to the unrestricted ( Table 2). The anisotropy of the relaxation rates in the restricted analysis is 2−3 times higher for the intramolecular couplings and 5 times higher for the couplings between water and peptide as compared to the corresponding anisotropies in the nonrestricted analysis ( Table 2).
The relaxation times T 1 and T 2 computed from the restricted analysis (Figure 7) are up to 78% shorter as compared to the unrestricted analysis (section 5.2). In the restricted analysis the anisotropic part makes up 40% of T 1 , while in the unrestricted analysis it is only 4%. The anisotropy in T 2 was considered negligible in the unrestricted analysis,  When the relaxation rates were computed by taking into account the water molecules only while they stay at the first hydration layer (Figure 6), a striking difference to the case of nonrestricted analysis (Figure 4) is seen: the anisotropies of the relaxation rates due to intramolecular couplings are not just different in magnitude but have different sign. This is likely due to the higher weighting given implicitly to the molecules which stay longer in the hydration layer.
The orientation dependence of T 1 (Figure 7) is inverted as compared to what was found for the unrestricted analysis (section 5.6) in the first hydration layer. This is due to the strong anisotropic contribution from the relaxation rate related to couplings between water and the peptide ( Figure 6) and the inverted anisotropy of the relaxation rates related to intramolecular couplings in water. However, not all of the TCFs of the water−peptide couplings were the shape of a multiexponential, which made the TCF fits less accurate than in the case of the analysis based only on the initial positions of the water molecules.
6.2. Comparison to Previous Experiments and Simulations. The correlation times of intramolecular couplings are directly the rotational correlation times of the H−H vector in a rigid water molecule. The 2.45 ps correlation times obtained for the bulk water correspond well to reorientation times in other simulation works at similar temperatures and pressures. 23,24,40 The correlation times of intramolecular couplings in the first hydration layer in the unrestricted analysis are roughly 3 times longer than the bulk water intramolecular correlation times. This translates to a rather moderate slowdown of water reorientation in the vicinity of the peptide. 23 In the restricted analysis, strikingly longer correlation times of intramolecular couplings, ranging from 39.3 to 150 ps were observed. However, this analysis gives more weight to those water molecules that stay longer times close to the peptide.
Earlier, by analyzing the residence times of water molecules in different hydrogen-bonding conditions, we found that most of the slowdown appears to result from cases in which the water molecules act as donors in a hydrogen bond. 41 The residence times and hydrogen-bonding has been analyzed for a similar collagen-like peptide more thoroughly earlier by Tourell and Momot. 17 Our T 1 = 3.42 s computed from the bulk water simulation is slightly shorter than the measured T 1 for pure water, which ranges from 3.57 s (25°C) 42 to T 1 = 3.838 s (27.50°C). 43 Previous molecular dynamics simulations using the same water model have given the value T 1 = 3.8 s at 298 K. 24 We were not able to trace the source of the difference to our result, although bulk simulations in different ensembles as well as removing the cutoff from the analysis of the bulk simulation were tested.
The correlation times computed from our simulations with the unrestricted analysis range from 2.45 to 36 ps (Table 1) so ω 0 τ c ≪ 1 and both R 1 and R 2 are to an excellent approximation independent of ω 0 , when ω 0 is below 10 9 1/s. Even in the case of the longest simulated correlation time from the restricted analysis, τ = 468 ps, the relaxation rates (and times) are independent of ω 0 = 2πν 0 , when ν 0 ≪ 340 MHz. Experimentally it has been found that T 1 in collagen solutions increases with the Larmor frequency of the proton, ν 0 , in the range ν 0 = 4.5, ···, 62 MHz. 44,45 This suggests that dipolar relaxation, at least as described with our simple model, is not the main source of longitudinal relaxation in articular cartilage. T 1 is also generally thought to be isotropic or weakly anisotropic in cartilage 4−7 and weakly anisotropic in tendon. 4,8 We found that T 1 was generally slightly anisotropic in the first hydration layer of the peptide but in practice this effect may be hard or impossible to observe in tissues such as articular cartilage, with an imperfect alignment of collagen fibers. When computation of the TCFs was strictly restricted to the molecules staying in the first hydration layer, the anisotropy of relaxation times was amplified and reversed due to the implicit stronger weighting of the molecules with longer residence times. In the case of the restricted analysis, the transverse relaxation time in the first hydration layer approaches those determined for the "bulk" water fraction in a solution with high (37 wt %) collagen concentration. 46 In our view, this water fraction would be better described as "fastexchanging", since all of the water has access to the collagen surface and forms short-lived hydrogen-bound states.
The transverse relaxation times in articular cartilage are typically in the range of tens of milliseconds when 400 2 0 = ω π MHz, 7 whereas even the shortest T 2 obtained in the restricted analysis of the first hydration layer in the current work was longer than 200 ms. If intramolecular couplings in water would be the dominating contributor to relaxation, according to eq 18, the experimentally observed relaxation times would correspond to correlation times in the subnanosecond (10 −10 s) range. Compared to the bulk water rotational correlation times, this would represent a slowdown of 2 orders of magnitude. In our work, only a moderate slowdown of the rotational correlation time was observed. T 2 is known to be strongly anisotropic in articular cartilage and has a maximum when the collagen fibrils are at 55°angle with respect to the main magnetic field. 6,10 Our simulated orientation-dependence of T 2 (Figure 7) and the oscillation factor ( Figure 5) suggest that the magic-angle effect at 55°c ould more likely be a dephasing effect caused by RDCs, than a result of slowdown of the reorientational motion of water. The former could be a result of anisotropic reorientation, whereas the latter is commonly attributed to water bound to the collagen as water bridges. 11, 12 6.3. Limitations and Outlook. We note that when the overall motion is anisotropic and residual dipolar interaction is included, the common formula [eq 18] for the transverse relaxation in liquids can only be taken as an approximation in anisotropic systems and a full treatment should include the dephasing effects of the RDC [eq 14]. Since the theory is strictly valid only for an isolated homonuclear spin-1/2 pair, the results for intermolecular couplings are estimates at best. Cross-correlations between different dipolar-coupled spin pairs were assumed negligible, but the assumption was not validated. A more severe limitation of the calculations was using a cutoff for the intermolecular couplings. This affects mostly the intermolecular RDCs. The even shorter cutoff used in the bulk simulation was not found to be a significant source of the error in the relaxation times of that system. Long-range corrections could be envisioned with, for example, Ewald summation or other methods for the RDCs. 47−50 Self-diffusion of water should be included with a realistic time scale, possibly unattainable in atomistic simulation. A larger system of bulk water would allow a longer cutoff.
Naturally, a single, relatively short triple-helix can not capture all aspects of the collagen−water interactions. For example, the water trapped inside fibrils could have significant interactions with several tropocollagen molecules at the same time. This shortcoming of the model could possibly be overcome by constructing a (staggered) bundle of parallel peptides. The amount of water inside the fibrils should be small as compared to the total water content in the system. The reorientation and diffusion of interstitial water could be significantly different from the relatively free water dynamically bonding with the collagen-like peptide. The shortness of the 1QSU peptide as compared to a real tropocollagen could overrepresent the interactions of water with the ends of the peptide. Moreover, the strongly charged termini in the peptide could also contribute to this overestimation. This aspect could be tested by, for example, excluding water near the termini from the analysis.
The current simulations are too short to realistically describe the movement of the peptide in solution, which is much slower than the movement of water molecules. Anyway, we are more interested in the cases in which the water molecules interact with a collagen strand that preserves its orientation, at least on average. The difference in the time scales of the water and peptide motions justifies separating them in the analysis.
Unfortunately, the simple peptide model can not accurately reproduce dephasing effects. A realistic description of water movement in an anisotropic collagen network would be needed to address these. Attempts in this direction have been carried out by others. 20 Altogether our work presents a first step toward modeling dipolar relaxation of water hydrogens in a guest−host system in which a collagen matrix acts as a host and water molecules as a guest. More realistic fibril models could be created by wrapping a tropocollagen model into a unit cell. 51,52 This would require a complete type-II tropocollagen molecule instead of the fragment-like peptide employed here. Approaching experimental conditions could be made by including effects of water diffusion inside and outside the fibrils as well as exchange with the bulk. Information on the orientation distribution of collagen fibrils in a tissue could be obtained, for example, via polarized light microscopy. 53

CONCLUSIONS
We found that in the first hydration layer of the peptide, the longitudinal relaxation time T 1 is weakly anisotropic with extrema found when the long axis of the peptide is oriented either perpendicular or parallel to the main magnetic field. The source of this anisotropy is mainly the intramolecular dipole− dipole couplings, with water−peptide couplings yielding weaker anisotropic contributions. Overall, T 2 was found to be isotropic, although individual contributions from the intramolecular couplings in water, as well as intermolecular couplings between 1 H in water and the peptide resulted in anisotropic relaxation rates. For water molecules further away from the peptide, both relaxation times are isotropic. When the water molecules were free to escape the hydration layers, the correlation times of the dipolar couplings were found to be mostly in the picosecond range. This makes the relaxation rates practically independent of B 0 . When only the molecules staying in the first hydration layer of the peptide were analyzed, the correlation times were significantly longer, but still well below one nanosecond. In this case of restricted analysis, the relaxation rates/times were found to be strongly anisotropic, with the exception of intermolecular couplings between water molecules. In light of our results the experimentally observed transverse relaxation rate in articular cartilage can not be explained by dipolar relaxation in the framework of the Redfield theory, but appears to be dominated by other relaxation mechanisms or dephasing effects resulting from RDCs.
Redfield relaxation theory for dipolar-coupled spin-1/2 pair with residual dipolar coupling, with cylindrically symmetric relative motion of the two nuclei; radial distribution functions between water and the peptide; justification of cylindrical symmetry; anisotropy of the simulated spectral densities (PDF)