Realtime optimization of multidimensional NMR spectroscopy on embedded sensing devices

The increasingly ubiquitous use of embedded devices calls for autonomous optimizations of sensor performance with meager computing resources. Due to the heavy computing needs, such optimization is rarely performed, and almost never carried out on-the-fly, resulting in a vast underutilization of deployed assets. Aiming at improving the measurement efficiency, we show an OED (Optimal Experimental Design) routine where quantities of interest of probable samples are partitioned into distinctive classes, with the corresponding sensor signals learned by supervised learning models. The trained models, digesting the compressed live data, are subsequently executed at the constrained device for continuous classification and optimization of measurements. We demonstrate the closed-loop method with multidimensional NMR (Nuclear Magnetic Resonance) relaxometry, an analytical technique seeing a substantial growth of field applications in recent years, on a wide range of complex fluids. The realtime portion of the procedure demands minimal computing load, and is ideally suited for instruments that are widely used in remote sensing and IoT networks.

S e e f T T dT dT where  is the experimental noise, θ is a calibration coefficient of the instrument, τ 1 and τ 2 are prescribed parameters that the spectrometer traverses through, with S(τ 1 , τ 2 ) the corresponding recorded signals. Inversion methods 22 are applied to obtain the sample {T 1 , T 2 } correlation spectrum, f(T 1 , T 2 ). As real-life samples often have broadly distributed f(T 1 , T 2 ), we made no assumptions other than T 1 ≥ T 2 23 on the mathematical constructs of the spectra.
Three classes of fluids were considered in the sequence design. As shown in Fig. 1, class A contained components that were longer than 0.1 s in both T 1 and T 2 dimensions; class B embodied high T 1 /T 2 ratios, where T 1 had components longer than 0.1 s while T 2 spanned [1 ms, 0.1 s]; and class C had relatively short T 1 and T 2 that each spanned [1 ms, 0.1 s]. Accordingly, Table 1 shows the optimal sequences for the respective fluid classes (i.e. sequence α for fluid A, β for B, and γ for C). In particular, both sequences α and β had τ 1 up to 10 s, capable of measuring T 1 up to 2 s, while sequence γ had the maximal τ 1 of 1 s that sufficed to measure T 1 up to 0.2 s. Meanwhile, sequence α had the maximal τ 2 = N e × t e of 10 s, capable of measuring T 2 up to 2 s, in contrast to sequences β and γ with the maximal τ 2 of 0.6 s. The shorter echo spacing, t e , used in sequences β and γ than in sequence α could further help resolve fast T 2 components (Fig. S1).  Table 1. The three optimal pulse sequences. PS is short for pulse sequence, t e is echo spacing, N e is the number of echoes, τ 1,min is the minimal τ 1 , τ 1,max is the maximal τ 1 , N 1 is the number of τ 1 , WT is the polarization time, and NA is the number of acquisitions. t 90 and t 180 are respectively 90° and 180° pulse lengths that are experimentally determined, and runtime is the total experimental time. www.nature.com/scientificreports www.nature.com/scientificreports/ Ideally, sequences shall always be applied to the intended samples under study; but in continuous measurements on samples of changing properties, occasions do arise in which a sequence is suboptimally applied. As sensors generally couldn't foresee temporal progression of sample properties, any combinations of fluid class and pulse sequence are practically probable. Here we show applications of each sequence to three exemplary fluids, namely dodecane (fluid A), emulsified fluid (fluid B), and glycerol (fluid C), in Fig. 2. In the 3 by 3 panels, the diagonal time-domain images were acquired by the respective optimal sequences, while the measurements that corresponded to the off-diagonal ones were either inefficient or erroneous (Fig. S1). The key challenge was to discern the fluid class from live time-domain images, regardless of the sequences in use, and apply the intended www.nature.com/scientificreports www.nature.com/scientificreports/ one in the subsequent runs. In practice, we used three trained ECOC (error-correcting output codes) learners 24 , a class of supervised learning models, for the realtime multiclass classification and inference task.
Supervised learning requires large quantities of labeled datasets for model training. The training sets can consist of either prior measurements on samples of known properties or forward-modeled simulations. We opted for the later approach thanks to the well-defined functional form of Eq. 1. Specifically, we approximated probable T 1 − T 2 distributions by a large ensemble of synthetic distributions, each consisting of three components of randomly generated   T T { , } 1 2 pairs. Since the measurement volume was a constant and filled with fluids of similar proton density, we assigned each component by a randomly-generated weighting coefficient, μ , that sums to unity. The time-domain data, S T , for model training were generated as: T e e n n n T N t T . After calibrating the experimental setup, we set θ to 1.85 and T  to a Gaussian noise with zero mean and 0.012 standard deviation.
To generate T T { , } 1 2   pairs, we stochastically sampled in a two-dimensional space, where each dimension consisted of 100 logarithmically distributed numbers from 1 ms to 2 s and all pairs satisfied the relation T T 1 2 ≥   . In total, 30,000 −   T T 1 2 distributions were created. Subsequently, we sifted the sampled distributions, one by one, through a set of classification criteria, and labeled them accordingly (Fig. S3). A given distribution was labeled class A if the longest  T 2 > 0.1 s and its associated weighting coefficient ≥0.05, labeled class B if the longest  T 2 < 0.1 s, its associated T 1  > 0.15 s, and its associated weighting coefficient ≥0.05, and labeled C if the longest  T 1 < 0.1 s. As a result, 11,663 were assigned to class A, 11,835 to B, and 1449 to C. www.nature.com/scientificreports www.nature.com/scientificreports/ To further reduce the size of training datasets, we exploited the separable structures of the functional form, and applied singular value decomposition (SVD) on T 1 and T 2 kernels independently 22 . Consequently, for any given sample, the size of compressed datasets was 1.57 KB when acquired by sequence α, 1.34 KB when acquired by β, and 1.25 KB when acquired by γ. As elaborated in the supplementary information, a near 1000-fold reduction in memory usage was achieved with the SVD compression.
Subsequently, we trained three ECOC classifiers, an ensemble method for multiclass classification problem 24 ; each classifier is used while running the corresponding pulse sequence. The ECOC models encoded the binary classification results from three linear support vector machines (SVM) into a coding design matrix 25 , using the "one-versus-one" strategy that distinguished a pair of labeled time-domain patterns in the training set while ignoring the third fluid class. For sequences α, β, and γ, the ECOC classifiers had the respective size of 4.7, 4.1 and 3.8 KB. In total, less than 13 KB of fast memory was required to store the models.
In the inference stage of classifying a new 2D NMR dataset, we utilized the loss-weighted decoding scheme 26 to aggregate predictions of the binary learners, in which the weighted "hinge" error functions 25 over all binary losses were minimized. After running a pulse sequence, the number of floating-point calculations for classifying the generated data, after normalization and SVD compression, is fewer than 700. More details of the model training, validation and inference can be found in the supplementary information.
We performed realtime optimizations of continuous NMR experiments with the trained ECOC classifiers, as illustrated in Fig. 3. The mobile NMR sensor 27 , shown in Fig. 3B, was miniaturized largely due to the use of an NMR ASIC (Application Specific Integrated Circuit) 3 . The NMR probe, embodied in a Halbach-array magnet, was made of a solenoid coil wound around a polymer capillary, interrogating fluid samples of 17 μL in volume. During operation, the spectrometer executed a selected pulse sequence that was downloaded from the laptop, on which the acquired data were input to ECOC classifiers for realtime inference. Figure 3C shows a series of experiments on sequentially displaced samples. The samples were six water/glycerol mixes of varying volume fractions and one emulsified fluid. As the volume fraction of glycerol increases, the relaxation times of the mixtures shorten from T 2 = T 1 = 2 s of pure water to T 2 = T 1 /2 = 15 ms of pure glycerol; meanwhile, the T 1 /T 2 ratio also inflates thanks to the elevated fluid viscosity. We started with sample 1 of pure water while applying sequence γ; the classifier correctly identified the fluid as type A and instructed to use sequence α for the subsequent run. Thereafter, sequence α was optimally applied for samples 1, 2 and 3. Sample 4 had T 1 /T 2 slightly above 1 with T 2 = 0.1 s. Consequently, the optimization routine signified the sample as type B. Sequence β was properly applied till sample 5 was loaded, which was classified as fluid C. Subsequently, sequence γ was applied on sample 5 and 6 of rather short relaxation times. Finally, we displaced glycerol by a well-gelled emulsion sample, which the ECOC classifier correctly deduced as a class B fluid; the spectrometer subsequently applied sequence β for the rest of the experiments.
In addition to physical displacement, a given sample could also evolve over time. For example, the emulsion fluids, which are multiphasic mixes of oil, brine, organoclays, and naturally-mined barite particles, could experience phase separation under static conditions. As the emulsion collapsed and solid particles sedimented, the emancipated oil exhibited a characteristic T 2 time much longer than the original fluid, calling for a different optimal sequence. Experimentally, we performed the optimization routine for continuously monitoring an emulsion sample under static conditions. As shown in Fig. 3D, the initial well-gelled emulsion fell in the class B fluid, to which sequence β was optimally applied. At about 50,000 s after the experiment commenced, signals of bulk oil started to appear with a T 2 at ca. 0.5 s. As volume fraction of free oil increased, the fluid gradually morphed to class A, with the corresponding optimal sequence α. Notably, A transition window presented at ca. 60,000 s, where the free fluid content was still marginal while the inference results hinged partially on noise realizations of each measurement.
In practice, it is important to ensure that each sequence has sensitivities over the entire numerical domain of T 1 − T 2 spectra under consideration. Failure to meet the requirement could cause misclassification and thereby erroneous results. For example, the maximum τ 1 in sequence γ, which is optimized for samples of fast relaxation times, should be designed so that the signal decays significantly with the maximal T 1 . Mathematically, it should satisfy 1 − exp(−τ 1,max /T 1,max ) ≫ σ, where σ 2 is the variance of the Gaussian noise. The relation is indeed held in the work, as τ 1,max = 1 s for sequence γ, T 1,max = 2 s, and σ 2 = 0.012 2 .
Although we focus on relaxometry, the method can be extended to other types of NMR measurements of increasing complexities, such as multidimensional spectroscopy and MRI, at the core of which are forward models of similar mathematical constructs (exponential, sine and cosine functions). In conjunction with minimal requirements on computing resources, the demonstrated approach may further NMR methods to a substantially broadened usage in a wide range of field applications.