Performance of DMEx models

The chemical exchange master equation is only homogeneous and analytically integrable in the simplest of cases and can acquire nonlinearities when one considers reversible exchange between distinguishable ensembles. The DMEx equations of motion presented here are fully converged in the exchange interaction and handle that interaction exactly but are not exact for the full evolution, as there is no preservation of the time ordering between the quantum and exchange degrees of freedom that would be present in a full solution. This could be approached by calculation of higher-order time derivatives but at an additional computational cost. However, the next section focuses on the difference between the traditional master equation approach (Eq. 7) and the DMEx solutions, which shows marked improvements in permissible time steps, and then compares the DMEx to QMC simulations where these time orderings are preserved, but the computational overhead is very large. Therefore, we will evaluate the performance of DMEx methods using the solution(t+t)=U(t)U+qtqexp(t2q)(KqE)(U(t)U)(30)

In this equation, Uexp(iH0t). This is an ideal computational method, as it only involves forward propagation of the solution, requires the fewest number of matrix operations, and produces linear evolution under the spin Hamiltonian and evolution to all orders in the exchange interaction. Equation 30 has a small intrinsic error associated with the solution, in that the first step only evolves quantum mechanically. A more accurate way to solve the equation of motion would be to evolve the initial density matrix backward in time by t/2 and then using Eq. 30 to generate the solution. Doing so shifts the actions of U and Kq by a half step and corrects for this initial error. However, we have found that this makes little difference in the solution; thus, we retain Eq. 30 so to avoid generating a nonintegral number of steps. To isolate errors arising from exchange, we have constructed all of the following simulations in the Hilbert space of the system, where one can exactly evaluate quantum evolution in systems up to 15 coupled spins.

Dynamic nuclear magnetic resonance spectra under pseudorotation have been studied and understood for decades. Spectral features are well resolved in the limit of slow exchange, which broaden and coalesce as the exchange rate increases, and ultimately result in line narrowing in the fast exchange limit. This is reflected in spectra of s-trioxane (9) undergoing ring inversion (Fig. 2A), where the axial (blue) and equatorial (red) have different chemical shifts and the geminal 2JHH coupling is observable. As exchange increases, the spectrum collapses to a singlet, as the axial and equatorial positions become, on average, magnetically equivalent. Similar effects appear for the tert-butyl rotation in t-BuPCl2 (Fig. 2B), which additionally exhibits a transition that is invariant under exchange and thus does not broaden (33).

We use s-trioxane ring inversion to model G2 pseudorotation (A) and tert-butyl rotation in t-BuPCl2 to model G3 pseudorotation (B) The graphs at the bottom compare the root mean square deviation (RMSD) of the generated magnetization as a function of the time evolution such as in Fig. 1C using t = 10 s (which is taken as the ground truth). The DMEx models retain good fidelity with no additional computational overhead, even with step sizes commonly being 10 times larger than were possible with traditional solutions.

For either of these systems, the pseudorotation matrices are generated by expressing a spin label permutation matrix in the appropriate basis. For convenience, we will use the Zeeman basis in this example. In the case of s-trioxane, where the axial and equatorial protons interchange, it is most convenient to setup the system such that axial protons have odd indices and equatorial protons have even indices. Then, the rotation R is given byR=P12P34P56(31)

Pij is the permutation matrix that interchanges the and states, whereas the and states are invariant under exchange. Using this method, it is trivial to arbitrarily reindex the entire system and is computationally efficient because the transformation from the original basis to the reindexed basis is unitary.

When calculating these spectra, we find that the traditional implementation and the DMEx converge to the same solution as t dt. However, the DMEx exhibits a substantially smaller error at any step size than the traditional implementation and only accrues an error on the order of 1% when the step size exceeds the average lifetime. In this limit, the traditional implementation loses stability and the trace of deviates from unity. This immediately provides the ability to take larger step sizes with the DMEx implementation. In the case of s-trioxane, an error in the solution of 1% requires t = 1 ms in the DMEx and t = 0.1 ms using the traditional implementation, thus requiring one to sample far fewer data points. In considering all moments of the exchange interaction, the radius of convergence of the Dyson expansion is far larger than it would be by assuming conditions similar to those used for exchange.

While these model systems provide illustrative examples of the performance of the DMEx model, they are far from the more challenging cases in dynamic systems. As noted previously, an interesting system that has gained much attention in the past decade is the hyperpolarization method SABRE, wherein large nonthermal nuclear magnetization is distilled from parahydrogen via reversible interactions with an iridium catalyst. Current efforts are focused on optimizing the extraction of spin order from parahydrogen, which requires accurate modeling of the quantum and exchange dynamics in realistic systems. For reference, an example simulation of the coupled coherent and exchange dynamics that drive SABRE hyperpolarization is shown in Fig. 3A, where the evolution of the 15N polarization is calculated under the experimental conditions for SABRESHield Enables Alignment Transfer to Heteronuclei (SHEATH) (16).

An example of SABRE hyperpolarization dynamics is shown for reference, calculated on a six-spin 15N SABRE-SHEATH system (A). Comparing the DMEx and a QMC treatment, which is viewed as the gold standard but is computationally inefficient, indicates that there is a genuine but small difference of 0.142 0.018%, on average, between the two solutions [red data, (B)]. The convergence error of the QMC is indicated by the black line. This error in the DMEx is attributed to the loss of the time orderings between the quantum and exchange degrees of freedom that are retained in the QMC. Even with nonlinear effects incorporated in the simulation, the DMExFR2 exhibits a larger radius of convergence over the traditional implementation by approximately a factor of 4 (C) and an improved self-consistency (parameterized by k, the error in the predicted exchange rate) (D), which uses t = 10 s solutions as ground truth.

In deriving the DMEx, we began with the ansatz that exchange could be considered as a time-continuous perturbation of the ensemble, but it is interesting to see when this assumption fails. The perturbation generated by exchange in the slow exchange limit is small, allowing the solution to be largely dictated by the quantum dynamics, and conversely, in the fast exchange limit, quantum evolution cannot generate large excursions from equilibrium when constantly disrupted by exchange. In the intermediate regime where SABRE exists, characterized by exchange rates on the order of the dominant couplings, it is no longer trivial to motivate that large excursions from equilibrium would not be impactful on the dynamics. To probe this, we compared our previous QMC model for SABRE (5) against the DMEx on a three-spin SABRE system (Fig. 3B) with a dominant coupling of 2JNH = 24 Hz. In this regime, there is a significant difference between the convergence error of the QMC solution (QMC) and the DMEx solution; however, this error is, on average, only 0.142 0.018%. Note that this analysis is limited to the smallest systems, given the large cost of iterating the QMC solution, and the error accrued by the DMEx is negligible on simulation time scale relevant to experimental guidance.

When modeling more complex systems, such as those often found in SABRE, it is critical for the cost of the DMEx to be augmented with an efficient method for exploring complex interactions, otherwise circumventing the benefits of an infinite-order treatment by excessive computational costs. In SABRE, these interactions include quantum evolution of multiple species, rebinding of previously polarized ligands to the activated complex, binding site competition with spin-inert coligands, and relaxation. We call this SABRE-specific model the DMExFR2 to indicate that free ligand, rebinding, and relaxation effects are included. The most efficient way of accomplishing an efficient implementation of exchange is by expressing the interactions as block diagonal with respect to individual manifolds, which we call manifold-diagonal for simplicity and will motivate using the example of SABRE.

In SABRE, we primarily consider two different species: one in which the hyperpolarization target is bound to the iridium, which we call bound species, and another in solution, which we call free species. Coherent evolution in the manifolds is established by separately propagating a bound species density matrix (bS) and the dissociated free species density matrix (fS) under their respective nuclear spin HamiltonianstfS=i[HfS,fS](32)tbS=i[HbS,bS](33)

We begin with the dynamics of the free species. Exchange facilitates association of free species to the catalyst (Ka,fSfS), acting on fS to remove free species from the manifold as they bind the complex, and allows for dissociation of bound ligand Kd,fSbS, adding species back to the manifold. Both of these exchange processes happen at the exchange rate of the ligand, kN, with an action scaled proportional to the ratio of theconcentration of the iridium complex to the free ligand ([Ir]/[S]) to account for the inherent trace normalization of density matrices. The association operator is then simplyKa,bSfS=[Ir][S]fS(34)

The dissociation operator then deposits an equivalent number of ligands from the bound species subsystem into the free species subsystem. For the case where both available binding sites in the iridium complex are exchanging with the target ligand, there are distinct subsets of the nuclear spins in the bound species (Sa and Sb), which may dissociate to join the free species with equal probability. We average these possibilities to generate the dissociation operator for the free species, remembering to apply the concentration scaling factor for exchange between manifoldsKd,fSbS=[Ir]2[S](Tr{H2Sa}bS+Tr{H2Sb}bS)(35)

Combining Eqs. 34 and 35 yields the equation of motion for the free species with exchange interactionstfS=i[HfS,fS]+kN[Ir][S](Kd,fSbSfS)(36)where kNkNexp(tkN/2) and will be used as a notation for the DMEx rate going forward.

The bound species has two exchange interactions: one for the simultaneous exchange of a ligand and the hydrides occurring at rate of kH and one for the exchange of target ligands at a rate of kN kH. We will formulate the exchange operator for the bound species as a single entity, Kex, which takes multiple manifolds as arguments. Hydride exchange is restricted to occur only during ligand exchanges as the complex form a tetrahydride intermediate to facilitate this reaction. In the case where both parahydrogen and ligand exchange occur concurrently, we exchange the portion t(ka,H/kN) of bS to reflect the new hydride population and new ligand population. This may be written as(ka,HkN)pH2fSTr{H2,S}bS(37)where pH2 is the density matrix for pure singlet parahydrogen and Tr{H2, S} returns the density matrix for the ligand that remains bound. In the case where the hydrides do not exchange but the target ligand does, another portion of the density matrix t((1ka,H)/kN) must be reformulated to reflect the newly exchanged ligand(1ka,HkN)Tr{S}bSfS(38)where Tr{S}bS is the density matrix for the subsystem of the remaining ligand and parahydrogen. This projection must be constructed carefully to ensure that the coherences are appropriately retained between the hydrides and remaining ligand. Note that while we are exchanging between the free and bound species subsystems, the scaling factor is not needed as the free species density matrix is, by definition, trace normalized. Therefore, one free ligand equivalent leaving the free species will look like one free ligand equivalent associating with the bound species. As this free ligand leaves the free species though, the appropriate reduction in the free species density matrix must be scaled by the concentration ratios. The full exchange operators can now be written as a combination of these two componentsK{S}(bS,fS)=(1ka,HkN)Tr{S}bSfS+(ka,HkN)pH2fSTr{H2,S}bS(39)

The two possible ligand exchanges from the two available binding sites, a and b, then average together to give the final exchange operator for the bound speciesKex=12(K{Sa}+K{Sb})(40)

Note that Eqs. 35 and 40 contains terms that are quadratic in the magnetization density, arising from the effects of rebinding ligands that have already interacted with the species. Hence, this is a second-order nonlinear partial differential equation, which must be solved simultaneously with the equation of motion for the free species to define the full evolution of the system.

Furthermore, these nonlinearities are amplified as the ratio increases. It now becomes possible to efficiently represent the impact of concurrent evolution of the J-coupling networks in the free and bound species of the target ligand. In addition, we can now model the effects that various solution compositions will have on the polarization dynamics, given that rebinding of previously polarized ligand will significantly affect the evolution of the bound species under the nuclear spin Hamiltonian.

Even with the incorporation of the nonlinear terms to the DMEx, the solution convergence is still far faster than that of the traditional implementation (Fig. 3C), and the two models still converge in the limit when t dt. One can obtain the same error in the DMExFR2 with a step size that is four times larger than the traditional implementation. While we have focused on the accuracy of the simulation, its precision in reproducing input parameters, such as the exchange rate, are just as important, particularly as these models are used to extract physical parameters from experimental data. Under this condition, it is not only critical that the simulation is stable but also efficient, as large portions of phase space have to be searched to perform an experimental fit. To characterize the precision of the simulation, we introduce the parameter k, which defines the relative shift in the predicted exchange rate in the simulation (Fig. 3D). Unexpectedly, there was essentially no exchange rate at which the traditional implementation provided a solution that was stably precise. In contrast to that, the DMEx model essentially perfectly reproduces the input exchange rate until k 300 s1, and when nonlinearities are introduced to the simulation, the maximum deviation from the input exchange rate is only k 0.5%.

As noted previously, guided in silico exploration of novel experimental methods that increase the hyperpolarization of SABRE is the focus of optimization efforts in the community. With the improved stability of the DMEx models, it is possible to explore realistic systems with complex coupling networks and reduce the calculation to an obtainable cost by using large simulation step sizes (t > 1 ms). The flexibility of this formulation to be expressed in either Hilbert or Liouville space additionally provides access to much larger spin systems than previously possible.

The case of the canonical bis-(15N-pyridine) SABRE-SHEATH system is particularly interesting, as it contains 14 strongly coupled spins in just the iridium-bound manifold with 22 total spins and is perhaps the most prevalent system in 15N SABRE. As the full system is far outside the scope of previous exchange models for SABRE, it has been traditionally acceptable to truncate the spin system to an approximate system, fully or partially removing ancillary 1H nuclei, with the largest approximation reported in literature using a single 1H per ligand (4). Even in this case, the dynamics of the truncated model diverge greatly from the actual system dynamics (Fig. 4A), the latter which can be explicitly calculated using the DMExFR2 model with either 2-ms (black) or 5-ms (red) step sizes with only minor deviation between the solutions. As a result, the truncated model optimizes to exchange rates that are false while retaining a deviation of 10% from the actual system dynamics when reoptimized to the erroneous rates (Fig. 4B). This means that any physical parameters extracted from experimental data by the model will be greatly confounded by the truncation errors inherent to the formulation.

(A) To reduce computational overhead, virtually all calculations to date remove ancillary spins from the system, such as artificially reducing the number of protons on the pyridines in the bis-(15N) SABRE complex from 10 to 2 (A). Doing so alters the hyperpolarization dynamics (blue) as compared to the explicit 14-spin calculation (black), which is stable using the DMExFR2 models for step sizes even up to 5 ms (red). (B) Fitting the 14-spin calculated time evolution with this smaller model produces incorrect values of the exchange rates. (C) Including all relevant exchange pathways when modeling SABRE systems is also crucial for predicting accurate exchange parameters. Here, we fit the experimental (15N,13C)-acetonitrile hyperpolarization dynamics to DMExFR2 models with (solid) and without (dashed) coligand exchange effects. When neglecting these exchange pathways, the predicted exchange rates differ from the correct values by 44 to 92%.

To emphasize the efficiency and flexibility of this framework, we used the DMExFR2 model to fit the coherent hyperpolarization dynamics of (15N-13C)-acetonitrile when exciting the sample with short (milisecond) pulses tuned to a field near the SABRE resonance condition, as described in our prior work (5, 32). Coherent evolution is then interrogated by varying the resonant pulse length, which encodes the dynamics in the final polarization detected. This is a multicomponent SABRE system containing 21 total spins and requires consideration of hyperpolarization-inactive coligand effects to accurately describe the dynamics. These effects allow for additional exchange pathways to influence the dynamics of the system. One of the most critical ramifications arising in allowing the hyperpolarizable ligand to exchange between positions on the complex and thus with which parahydrogen-derived hydride the ligand is coupled. In the limit of fast exchange, this makes the hydrides appear equivalent and would prevent the singlet order from being converted into observable magnetization. When coligand effects are included (solid lines), the experimental data can be reproduced with high fidelity to experiment at multiple field conditions (Fig. 4C), such as when the resonant pulse is B = 1.65 (red) or 0.91 T (blue). Furthermore, the extracted exchange rates for these datasets are kN = 14.5 1.8 s1 and kH = 6.00 0.75 s1 for the B = 1.65 T data and kN = 15.0 3.3 s1 and kH = 4.50 0.98 s1 for the B = 0.91 T data. When coligand effects are neglected, the predicted exchange rates can have errors of 44 to 92%.

Properly simulating this system requires two seven-spin manifolds for the two conformers of the iridium complex, a five-spin manifold for the free (15N-13C)-acetonitrile, and a two-spin manifold for parahydrogen. Fitting the experimental data would be intractable within any of the previous formalisms for SABRE dynamics (3, 25) as a function of the system size. However, when built using the DMExFR2 model in conjunction with the manifold-diagonal formalism for exchange introduced here, each simulation dataset, which consists of 32 simulations lasting 30 s using t = 2 ms, requires approximately 15 min to calculate, making a grid optimization fit possible within a day.

More:

Infinite-order perturbative treatment for quantum evolution with exchange - Science Advances