On-the-fly, Sample-tailored Optimisation of NMR Experiments

Jonathan R. J. Yong and Mohammadali Foroozandeh*

Chemistry Research Laboratory, Department of Chemistry, University of Oxford,
12 Mansfield Road, Oxford, OX1 3TA, U.K.

*<EMAIL_ADDRESS>
NMR parameters must often be adjusted for maximum performance on different
samples or instruments. We show that, using our implementation of POISE
(Parameter Optimisation by Iterative Spectral Evaluation), a wide range of
experiments may be optimised on the fly in a fully automated fashion, using
only a few spectrum acquisitions.

###### Abstract

NMR experiments, indispensable to chemists in many areas of research, are
often run with generic, unoptimised experimental parameters. This approach
makes robust and automated acquisition on different samples and instruments
extremely challenging. Here, we introduce NMR-POISE (Parameter Optimisation by
Iterative Spectral Evaluation), the first demonstration of on-the-fly, sample-
tailored, and fully automated optimisation of a wide range of NMR experiments.
We illustrate how POISE maximises spectral sensitivity and quality with a
diverse set of 1D and 2D examples, ranging from HSQC and NOESY experiments to
ultrafast and pure shift techniques. Our Python implementation of POISE has an
interface integrated into Bruker’s TopSpin software, one of the most widely
used platforms for NMR acquisition and automation, allowing NMR optimisations
to be run without direct user supervision. We predict that POISE will find
widespread usage in academia and industry, where sample-specific and automated
experiment optimisation is mandatory.

Keywords: NMR spectroscopy $\cdot$ Analytical methods $\cdot$ NMR-POISE
$\cdot$ Optimisation $\cdot$ Automation

Nuclear magnetic resonance (NMR) spectroscopy is firmly established as one of
the most important analytical methods in chemistry[1, 2]. Although many NMR
experiments are designed to be as general as possible, the sheer diversity of
molecular structure means that experimental parameters must often be varied in
order to achieve optimal performance for a particular sample. Using poorly set
parameters can lead to inferior spectra, as measured by (for example) the
attainable signal-to-noise ratio (SNR) per unit time, or the presence of
artefacts arising from unwanted coherences. A typical way around this is to
use a “compromise” value, chosen to provide acceptable (albeit suboptimal)
performance across a range of samples.

In contrast, the direct optimisation of experimental parameters on a per-
sample basis is rarely carried out. When it is, it typically relies on trial-
and-error or otherwise inefficient methods, leading to highly tedious and
time-consuming processes. For example, the current best routines (e.g. popt in
Bruker’s TopSpin software) rely on systematic evaluation of parameter values
using an exhaustive grid search. A more desirable approach would be to utilise
an optimisation algorithm[3] which iteratively monitors the quality of a
series of spectra acquired with different parameter values, eventually
converging to optimal values for the sample under study. Although there exist
isolated examples of such iterative optimisations in laser[4], electron spin
resonance[5], and nuclear quadrupole resonance[6] spectroscopies, such a
procedure has only been used in NMR by Emsley and co-workers in the design of
dipolar decoupling pulses[7, 8, 9]. Note that the proposed method is distinct
from other model-based “optimisations”, e.g. for sparse sampling in
multidimensional NMR[10, 11] or relaxation measurements[12], where data from
multiple acquisitions are aggregated until a certain threshold, such as a
sufficient confidence in peak locations or extracted parameters, is reached.

Despite its apparent simplicity, iterative optimisation of NMR experiments
clearly enjoys little popularity. This can perhaps be attributed to a
perception that such a process would require too much time, as well as the
challenge of implementing it in a manner that can be suitably generalised for
multiple disparate contexts.

In this work, we dispel these notions by introducing NMR-POISE (Parameter
Optimisation by Iterative Spectral Evaluation), a method which enables the
efficient optimisation of a wide range of NMR experimental parameters and
yields demonstrably better spectra in only a matter of minutes. NMR-POISE is
provided as a free Python package, with a powerful user interface allowing
optimisations to be run directly from the command line of Bruker’s TopSpin
software. This means that optimisations can proceed in a fully automated
fashion without any user input, a benefit which is especially relevant in
high-throughput settings or when non-expert users are present. The
optimisation process (Figure 1) requires a cost function, a user-defined
function which measures how good a spectrum is; POISE provides a number of
built-in cost functions for typical use cases, e.g. for maximising or
minimising signal intensity. Currently, three optimisation algorithms are
provided: Py-BOBYQA, a trust-region interpolation method[13] (enabled by
default); the Nelder–Mead (NM) simplex method[14]; and the multidimensional
search (MDS) method[15].

Figure 1: Flowchart depicting the main steps in a POISE optimisation.

Experimental section

All NMR spectra in this work were acquired on Bruker AV III or NEO
spectrometers, with field strengths ranging from 400 to
$700\text{\,}\mathrm{MHz}$. POISE was executed on the spectrometer terminals
using Python 3.8.1, 3.8.5, and 3.9.2 on Windows 10. In general, POISE can be
run on any Bruker spectrometer; the minimum version of Python needed is 3.6.
Experimental details of individual optimisations are given in the Supporting
Information.

Results and discussion

The first example of a POISE control optimisation is the calibration of a
$None$ ^1H pulse, a quantity which is fundamental to practically every NMR
experiment[16]. In a simple ^1H pulse-acquire experiment, a perfect $None$
pulse would generate no detectable magnetisation. The intensity of the
resulting signal can thus be used as a cost function which we seek to
minimise.

Figure 2:  Optimisation trajectory for pulse width calibration on a sample of
$45\text{\,}\mathrm{m\textsc{m}}$ ferulic-acid in DMSO-$d_{6}$ using: 2
TopSpin’s popt, which uses a grid search; and 2 POISE, using the BOBYQA
algorithm. The grey line indicates the “true” objective function. Each
coloured arrow indicates a single function evaluation, unfilled dots represent
the initial point, and filled dots indicate optima at which the cost function
is minimised.

The benefits of direct optimisation, as compared to the grid search of popt,
can be easily appreciated upon consideration of Figure 2. By manually
acquiring an array of spectra with different pulse widths, we first
ascertained that the spectral intensity was minimised at
$48.4\text{\,}\mathrm{\SIUnitSymbolMicro s}$ (Figure S5); knowing this allows
us to evaluate the accuracy of the optimisation routines. popt faithfully
reproduces this optimum, but must evaluate a total of 41 spectra to search the
entire parameter space, which takes 5 minutes (Figure 2). Although the search
may in principle be terminated immediately after locating the minimum, this
requires constant monitoring by the user, since the optimum point is not known
a priori. On the other hand, POISE quickly homes in on the minimum without any
user intervention: for example, using the BOBYQA algorithm, it requires only
70 seconds and 9 spectrum acquisitions to identify an accurate minimum of
$48.3\text{\,}\mathrm{\SIUnitSymbolMicro s}$, even when given a suboptimal
starting point of $43\text{\,}\mathrm{\SIUnitSymbolMicro s}$ (Figure 2).
Convergence is even faster when a good initial point is provided (Tables S2,
S3 and S4). As a “rapid” alternative for pulse width calibration, the Bruker
pulsecal routine (which uses a nutation experiment[17]) completes in 37
seconds; however, it underestimates the $None$ null by almost
$2\text{\,}\mathrm{\SIUnitSymbolMicro s}$ (Table S2). For routine 1D ^1H NMR
spectra, the exact $None$ pulse width is often not of interest. However, by
instead searching for a maximum in the spectral intensity, POISE can likewise
be used to identify the pulse width corresponding to the Ernst angle[18]
(Section S5), which yields the maximal sensitivity for a specific repetition
time.

As long as a suitable cost function can be defined, POISE is capable of
providing rapid results with no compromise in accuracy. Our POISE interface is
thus designed to grant users great flexibility in defining their own cost
functions, which allows the optimisation of a wide range of parameters which
simply cannot be carried out using currently available tools. One such example
is the mixing time in nuclear Overhauser effect spectroscopy (NOESY)
experiments, which typically varies from hundreds of milliseconds to several
seconds. By finding the optimal mixing time for a given sample, the
sensitivity of the resulting NOESY spectrum can be maximised. Here, we perform
this optimisation using a selective 1D NOESY experiment (Figure 3), using the
intensity of the spectrum excluding the selectively irradiated peak as the
cost function. Starting from an initial value of $0.5\text{\,}\mathrm{s}$, a
typical value for many organic molecules, POISE is capable of finding the
optimum mixing time for the sample under study ($3.5\text{\,}\mathrm{s}$).
This optimisation converges in around 5 minutes (Table S8), leading to an
approximate twofold increase in NOESY crosspeak intensity for this particular
case (Figure 4).

Figure 3:  Pulse sequences used for POISE optimisations; 3 1D selective NOESY,
3 ASAP-HSQC,[19] 3 echo-planar spectroscopic imaging (EPSI) acquisition scheme
used in ultrafast NMR[20, 21], and 3 spin echo sequence used for evaluation of
PSYCHE[22] refocusing element, 3 1D NOESY sequence with presaturation, used
for solvent suppression.[23] Details of the optimisation processes, along with
complete lists of experimental parameters used, can be found in the Supporting
Information.

Where possible, POISE optimisations should be performed using 1D pulse
sequences, in order to minimise the optimisation duration. The NOE mixing time
optimisation shown here, for example, takes only a fraction of the time needed
for a 2D NOESY. With more dilute samples, a larger number of transients per
function evaluation may be required in order for a reliable cost function to
be measured; consequently, the optimisation as a whole will take a longer time
to run. However, since the same can be said of the final 2D spectrum that is
acquired, the benefits that can be reaped are correspondingly larger, making
the “return on investment” the same. A similar strategy can be used to search
for the “null” in an inversion-recovery (IR) experiment: by running the
optimisation on a 1D IR (i.e. a single increment of the full experiment),
$T_{1}$ relaxation times can be measured in a far shorter time than the full
2D IR (Section S7).

Figure 4:  2D NOESY spectra of a sample of $120\text{\,}\mathrm{m\textsc{m}}$
3-fluorophenylboronic acid in DMSO-$d_{6}$, obtained 4 before and 4 after
optimising the mixing time. Optimisation leads to an approximate $2\times$
improvement in crosspeak sensitivities. In order to accomplish this
sensitivity improvement, one would need to acquire the unoptimized experiment
with four times as many transients.

An example of a 2D experiment which can be optimised in a reasonable time is
the ASAP-HSQC experiment (Figure 3)[19], where the recovery delay between
transients is dramatically shortened. Decreasing the duration of the INEPT
delay leads to greater sensitivity through Ernst angle-style excitation of
^13C-bound proton magnetisation. Using the projection of the full 2D spectrum
onto the $f_{2}$ dimension as a cost function, POISE is able to optimise this
delay to approximately 60% of its original length. This takes approximately 3
minutes for the same sample as above (Table S12) and provides an 18–27%
increase in sensitivity (Figure 5). This procedure can be extended to the
optimisation of other delays in 2D spectra, particularly with “fast” variants
such as the ASAP method used here.[24]

Figure 5:  5 2D ASAP-HSQC spectrum of a sample of
$120\text{\,}\mathrm{m\textsc{m}}$ 3-fluorophenylboronic acid in DMSO-$d_{6}$.
5 $f_{2}$ projections of ASAP-HSQC spectra before (grey dotted line) and after
POISE optimisation of the INEPT delay (green solid line). The sensitivity
increases attained are indicated above each peak.

The logical extension of the preceding example is ultrafast 2D NMR
experiments, in which the combination of spatially encoded $t_{1}$ evolution
and the echo-planar spectroscopic imaging (EPSI) detection method enables
single-scan acquisition of 2D spectra.[20, 21] Despite the substantial efforts
of Giraudeau and coworkers in the development of an accessible protocol,[21]
obtaining high-quality ultrafast spectra remains a nontrivial task.

We illustrate here how POISE can be used to aid in the setup of ultrafast
experiments, in this case by compensating for systematic imbalances in the
positive and negative gradients used during EPSI acquisition. Such an
imbalance leads to a progressive shifting of echoes in the spatial frequency
($k$) domain (Figure 6), ultimately causing degradation of lineshapes in the
resulting spectra.[20] The pulse programme used for optimisation consists of
$None$ excitation, immediately followed by EPSI acquisition (Figure 3); the
amplitudes of the negative gradients are multiplied by a factor $\alpha$ in
order to compensate for any imbalance. Here, the cost function is the slope
generated by the echoes in a $(k,t_{2})$ data matrix (Figures 6 and 6), which
directly reflects the extent of gradient imbalance. On our system, POISE finds
the optimal value of $\alpha$ to be $1.00038$ in under a minute (Table S13).
Such a small imbalance may appear inconsequential, but in fact leads to marked
improvements in the resulting lineshapes of ultrafast spectra, such as the
TOCSY spectra shown in Figures 6 and 6.

Figure 6:  6 and 6: Magnitude-mode $(k,t_{2})$ data matrix obtained from the
excitation–EPSI pulse sequence (Figure 3), before and after optimisation of
the gradient imbalance factor $\alpha$. Data points acquired during the
negative gradients have been discarded. 6 and 6: The corresponding 2D
ultrafast TOCSY spectra of a sample of $45\text{\,}\mathrm{m\textsc{m}}$
ferulic acid in DMSO-$d_{6}$ ($10\text{\,}\mathrm{ms}$ MLEV-16 mixing), before
and after optimisation.

Finally, another major advantage of POISE is that its algorithms are capable
of simultaneously optimising multiple parameters in an efficient manner.
Although this is possible with a grid search, the time required scales
exponentially with the number of parameters, quickly rendering such a
procedure unfeasible. We show how this capability of POISE can be used to
search for optimal performance in PSYCHE broadband homonuclear decoupling[22].
Here, homodecoupling is accomplished by the application of a J-refocusing
element (JRE), comprised of two saltire chirped pulses[25] and a weak
gradient. These can be expressed in terms of several parameters, namely the
flip angle $\beta$, gradient amplitude $g$, bandwidth $\Delta F$, and duration
$\tau_{\mathrm{p}}$.

We used a 1D experiment for optimisation of these parameters, which consists
of a $None$ excitation pulse followed by a spin echo with the JRE in the
middle (Figure 3). Any non-idealities in the JRE can be captured in a
“spectral difference” cost function (Section S10), which compares this against
a standard pulse-acquire experiment. The simultaneous optimisation of all four
parameters can be completed in around 15 minutes and leads to noticeably
better decoupling performance for the sample under study, especially for the
peak at $1.36\text{\,}\mathrm{ppm}$ (Figure 7). Furthermore, the optimised
spectrum exhibits improved suppression of the artefacts at
$2.5\text{\,}\mathrm{ppm}$ (Figure 7), which arise from a challenging case of
strong coupling. Although the optimised spectrum has a decreased SNR, the cost
function used here also penalises sensitivity losses, thereby ensuring that
any SNR loss is not excessive.

Figure 7:  Insets from spectra of a sample of
$45\text{\,}\mathrm{m\textsc{m}}$ andrographolide in DMSO-$d_{6}$; artefacts
in PSYCHE spectra are indicated with asterisks. Left column: Ordinary coupled
^1H spectrum; middle column: unoptimised PSYCHE spectrum; right column: PSYCHE
spectrum after simultaneous optimisation of four parameters. 7 Peak at
$1.36\text{\,}\mathrm{ppm}$, where recoupling artefacts are eliminated through
optimisation. 7 Peak at $2.5\text{\,}\mathrm{ppm}$, where optimisation leads
to a substantial reduction in strong coupling artefacts.

In a similar fashion, we applied POISE to the 1D NOESY / presaturation
experiment commonly used for water suppression (Figure 3)[23], where up to
four parameters (transmitter offset, presaturation power, mixing time, and
presaturation duration) are optimised. This procedure takes between 4 and 25
minutes, depending on the number of parameters under optimization (Section
S11). By applying appropriate upper bounds on the presaturation power and
duration, POISE greatly reduces the intensity of the residual water peak,
while also ensuring that nearby peaks of interest are not lost (Figure 8).
Finally, a further example of multiple-parameter optimisation is presented in
the context of diffusion NMR, enabling a fully automated diffusion ordered
spectroscopy (DOSY) experiment (Section S12).

Figure 8:  Spectra of a sample of rodent urine. Left column: Unoptimised
spectrum; middle column: spectrum following optimisation of two parameters
(transmitter offset and presaturation power); right column: spectrum after
optimisation of all four parameters. 8 Peak at $5.15\text{\,}\mathrm{ppm}$. 8
Residual water peak. 8 Peak at $4.56\text{\,}\mathrm{ppm}$. Note that the
flanking peaks in 8 and 8 are largely unaffected by the improvement in water
suppression.

Through a series of examples, we have shown how POISE can be used to great
effect in the optimisation of NMR parameters. Each optimisation provides
parameter values that are tailored to the instrument and sample under study,
typically in only a few minutes and in a fully automated fashion through the
TopSpin command line or acquisition scripts. Importantly, POISE allows users
to define arbitrary cost functions, making it a highly versatile tool that can
be adapted to the myriad uses of NMR in modern chemistry. POISE is fully
compatible with, and will be immediately usable in, other fields such as
solid-state and biomolecular NMR. More generally, we also note that the
iterative optimisation process demonstrated here can be extended to other
forms of spectroscopy.

The source code of NMR-POISE is publicly hosted on GitHub:
https://github.com/foroozandehgroup/nmrpoise. Detailed instructions on
installation, cost function definition, and optimisation setup may be found in
the user documentation, which is included in the Supporting Information and
can also be found online at https://foroozandehgroup.github.io/nmrpoise.

## Acknowledgements

We thank Coralia Cartis, Tim Claridge, Jean-Nicolas Dumez, David Goodwin, Iain
Swan, and Chris Waudby for helpful discussions; as well as Fay Probert and Abi
Yates for providing the rodent urine sample. J.R.J.Y. thanks the Clarendon
Fund (University of Oxford) and the EPSRC Centre for Doctoral Training in
Synthesis for Biology and Medicine (EP/L015838/1) for a studentship,
generously supported by AstraZeneca, Diamond Light Source, Defence Science and
Technology Laboratory, Evotec, GlaxoSmithKline, Janssen, Novartis, Pfizer,
Syngenta, Takeda, UCB, and Vertex. M.F. thanks the Royal Society for a
University Research Fellowship and a University Research Fellow Enhancement
Award (Grant Nos. URF\R1\180233 and RGF\EA\181018).

## References

  * [1] John Cavanagh  “Protein NMR Spectroscopy: Principles and Practice”  Burlington, Mass.: Academic Press, 2007
  * [2] Timothy D.. Claridge  “High-Resolution NMR Techniques in Organic Chemistry”  Amsterdam: Elsevier, 2016
  * [3] Jorge Nocedal and Stephen J. Wright  “Numerical Optimization”  New York: Springer, 2006
  * [4] Christopher J Bardeen et al.  “Feedback quantum control of molecular electronic population transfer”  In _Chem. Phys. Lett._ 280, 1997, pp. 151–158  DOI: 10.1016/S0009-2614(97)01081-6
  * [5] David L. Goodwin, William K. Myers, Christiane R. Timmel and Ilya Kuprov  “Feedback control optimisation of ESR experiments”  In _J. Magn. Reson._ 297, 2018, pp. 9–16  DOI: 10.1016/j.jmr.2018.09.009
  * [6] Cristian Monea  “Optimization of NQR excitation sequences using black-box techniques”  In _J. Magn. Reson._ 321, 2020, pp. 106858  DOI: 10.1016/j.jmr.2020.106858
  * [7] Gaël De Paëpe, Paul Hodgkinson and Lyndon Emsley  “Improved heteronuclear decoupling schemes for solid-state magic angle spinning NMR by direct spectral optimization”  In _Chem. Phys. Lett._ 376, 2003, pp. 259–267  DOI: 10.1016/S0009-2614(03)00966-7
  * [8] Bénédicte Elena, Gaël Paëpe and Lyndon Emsley  “Direct spectral optimisation of proton–proton homonuclear dipolar decoupling in solid-state NMR”  In _Chem. Phys. Lett._ 398, 2004, pp. 532–538  DOI: 10.1016/j.cplett.2004.09.122
  * [9] Elodie Salager et al.  “Homonuclear dipolar decoupling with very large scaling factors for high-resolution ultrafast magic angle spinning ^1H solid-state NMR spectroscopy”  In _Chem. Phys. Lett._ 498, 2010, pp. 214–220  DOI: 10.1016/j.cplett.2010.08.038
  * [10] Hamid R. Eghbalnia et al.  “High-Resolution Iterative Frequency Identification for NMR as a General Strategy for Multidimensional Data Collection”  In _J. Am. Chem. Soc._ 127, 2005, pp. 12528–12536  DOI: 10.1021/ja052120i
  * [11] Alexandar L. Hansen and Rafael Brüschweiler  “Absolute Minimal Sampling in High-Dimensional NMR Spectroscopy”  In _Angew. Chem. Int. Ed._ 55, 2016, pp. 14169–14172  DOI: 10.1002/anie.201608048
  * [12] Y.-Q. Song, Yiqiao Tang, M.. Hürlimann and D.. Cory  “Real-time optimization of nuclear magnetic resonance experiments”  In _J. Magn. Reson._ 289, 2018, pp. 72–78  DOI: 10.1016/j.jmr.2018.02.009
  * [13] Coralia Cartis, Jan Fiala, Benjamin Marteau and Lindon Roberts  “Improving the Flexibility and Robustness of Model-based Derivative-free Optimization Solvers”  In _ACM Trans. Math. Softw._ 45, 2019, pp. 1–41  DOI: 10.1145/3338517
  * [14] J.. Nelder and R. Mead  “A Simplex Method for Function Minimization”  In _The Computer Journal_ 7, 1965, pp. 308–313  DOI: 10.1093/comjnl/7.4.308
  * [15] J.. Dennis and Virginia Torczon  “Direct Search Methods on Parallel Machines”  In _SIAM J. Optim._ 1, 1991, pp. 448–474  DOI: 10.1137/0801027
  * [16] Paul A. Keifer  “$90^{\circ}$ pulse width calibrations: How to read a pulse width array”  In _Concepts Magn. Reson._ 11, 1999, pp. 165–180  DOI: 10.1002/(SICI)1099-0534(1999)11:3¡165::AID-CMR4¿3.0.CO;2-D
  * [17] Peter S.. Wu and Gottfried Otting  “Rapid pulse length determination in high-resolution NMR”  In _J. Magn. Reson._ 176, 2005, pp. 115–119  DOI: 10.1016/j.jmr.2005.05.018
  * [18] R.. Ernst and W.. Anderson  “Application of Fourier Transform Spectroscopy to Magnetic Resonance”  In _Rev. Sci. Instrum._ 37.1, 1966, pp. 93–102  DOI: 10.1063/1.1719961
  * [19] David Schulze-Sünninghausen, Johanna Becker and Burkhard Luy  “Rapid Heteronuclear Single Quantum Correlation NMR Spectra at Natural Abundance”  In _J. Am. Chem. Soc._ 136, 2014, pp. 1242–1245  DOI: 10.1021/ja411588d
  * [20] Lucio Frydman, Adonis Lupulescu and Tali Scherf  “Principles and Features of Single-Scan Two-Dimensional NMR Spectroscopy”  In _J. Am. Chem. Soc._ 125, 2003, pp. 9204–9217  DOI: 10.1021/ja030055b
  * [21] Boris Gouilleux, Laetitia Rouger and Patrick Giraudeau  “Ultrafast 2D NMR: Methods and Applications”  In _Annu. Rep. NMR Spectrosc._ , 2018, pp. 75–144  DOI: 10.1016/bs.arnmr.2017.08.003
  * [22] Mohammadali Foroozandeh et al.  “Ultrahigh-Resolution NMR Spectroscopy”  In _Angew. Chem. Int. Ed._ 53, 2014, pp. 6990–6992  DOI: 10.1002/anie.201404111
  * [23] Ryan T. Mckay  “How the 1D-NOESY suppresses solvent signal in metabonomics NMR spectroscopy: An examination of the pulse sequence components and evolution”  In _Concepts Magn. Reson._ 38A.5, 2011, pp. 197–220  DOI: 10.1002/cmr.a.20223
  * [24] Paul Schanda  “Fast-pulsing longitudinal relaxation optimized techniques: Enriching the toolbox of fast biomolecular NMR spectroscopy”  In _Prog. Nucl. Magn. Reson. Spectrosc._ 55, 2009, pp. 238–265  DOI: 10.1016/j.pnmrs.2009.05.002
  * [25] Mohammadali Foroozandeh  “Spin dynamics during chirped pulses: applications to homonuclear decoupling and broadband excitation”  In _J. Magn. Reson._ 318, 2020, pp. 106768  DOI: 10.1016/j.jmr.2020.106768

Supporting Information

On-the-fly, Sample-tailored Optimization of NMR Experiments

Jonathan R. J. Yong and Mohammadali Foroozandeh*

Chemistry Research Laboratory, Department of Chemistry, University of Oxford,
12 Mansfield Road, Oxford, OX1 3TA, U.K.

*<EMAIL_ADDRESS>
###### Contents

  1. S1 Raw data and reproducibility
    1. S1.1 Datasets
    2. S1.2 Figures
  2. S2 Routines
    1. S2.1 poise_1d AU programme
    2. S2.2 poise_1d_noapk AU programme
    3. S2.3 poise_2d AU programme
  3. S3 1D ^1H spectra for all samples
  4. S4 $\mathbf{90^{\circ}}$ pulse width calibration
    1. S4.1 Cost function
    2. S4.2 Reference grid search
    3. S4.3 Optimisation
    4. S4.4 Using POISE under automation
      1. S4.4.1 AU programme
      2. S4.4.2 Python script
  5. S5 Ernst angle
    1. S5.1 Cost function
    2. S5.2 Theoretical optima via $T_{1}$
    3. S5.3 Optimisation
    4. S5.4 Variable flip angle pulse-acquire pulse programme
  6. S6 NOE mixing time
    1. S6.1 Cost function
    2. S6.2 Reference grid search
    3. S6.3 Optimisation
    4. S6.4 1D NOE pulse programme
  7. S7 Inversion-recovery
    1. S7.1 Cost function
    2. S7.2 2D inversion-recovery
    3. S7.3 Optimisation
    4. S7.4 1D inversion-recovery pulse programme
  8. S8 ASAP-HSQC
    1. S8.1 Cost function
    2. S8.2 Reference grid search
    3. S8.3 Optimisation
    4. S8.4 ASAP-HSQC pulse programme
  9. S9 EPSI
    1. S9.1 Cost function
    2. S9.2 Reference grid search
    3. S9.3 Optimisation
    4. S9.4 Excitation–EPSI pulse programme
  10. S10 PSYCHE
    1. S10.1 Cost function
    2. S10.2 Flip angle reference grid search
    3. S10.3 Optimisations
    4. S10.4 AU programme used for optimisations
    5. S10.5 Python script for generating double saltire
    6. S10.6 1D spin echo pulse programme
  11. S11 1D NOESY-based solvent suppression
    1. S11.1 Cost function
    2. S11.2 Transmitter offset reference grid search
    3. S11.3 Optimisations
    4. S11.4 1D NOESY pulse programme
  12. S12 Diffusion NMR
    1. S12.1 Sequential optimisation strategy
    2. S12.2 Python script for automation
    3. S12.3 Sequential optimisation details
    4. S12.4 Simultaneous optimisation strategy
    5. S12.5 Oneshot DOSY pulse programmes
      1. S12.5.1 2D Oneshot
      2. S12.5.2 1D Oneshot
  13. S13 Documentation

Note: The Supporting Information sections mainly contain summaries of the
results presented in this paper, as well as technical details for reproducing
these results. It does not explain how to set up POISE, or to how to create
your own optimisation routines. For step-by-step instructions on these, please
refer to the documentation. A copy of it is attached at the end of the
Supporting Information, but the most updated version will be accessible at
https://foroozandehgroup.github.io/nmrpoise (a PDF copy is located at
https://foroozandehgroup.github.io/nmrpoise/poise.pdf).

##  S1 Raw data and reproducibility

All raw data used for this paper, as well as the Python scripts to generate
the plots in the main paper and SI, can be downloaded from Zenodo at the DOI
10.5281/zenodo.4698423.

###  S1.1 Datasets

The datasets that were used for each section are tabulated in Table S1.

Section | Description | Folder name | Expnos
---|---|---|---
$None$ pulse width | Reference grid search | P1_Scan | 1–161
| Optimisation (initial value $48\text{\,}\mathrm{\SIUnitSymbolMicro s}$) | P1_Opt | 4
| Optimisation (initial value $43\text{\,}\mathrm{\SIUnitSymbolMicro s}$) | P1_Opt | 5
| Optimisation (initial value $53\text{\,}\mathrm{\SIUnitSymbolMicro s}$) | P1_Opt | 6
Ernst angle | 2D inversion-recovery | Invrec_Opt | 2
| Optimisation (6 to $8\text{\,}\mathrm{ppm}$) | Ernst_Opt | 1001
| Optimisation ($6.79\text{\,}\mathrm{ppm}$ only) | Ernst_Opt | 1002
NOE mixing time | Reference grid search | NOESY_OptScan | 101–159
| Optimisation | NOESY_OptScan | 3
Inversion-recovery | 2D inversion-recovery | Invrec_Opt | 2
| Optimisation (6 to $8\text{\,}\mathrm{ppm}$) | Invrec_Opt | 3
| Optimisation ($7.08\text{\,}\mathrm{ppm}$ only) | Invrec_Opt | 4
ASAP-HSQC | Reference grid search | HSQC_OptScan | 101–133
| Optimisation | HSQC_OptScan | 2
EPSI | Reference grid search | EPSI_Scan | 95–116
| Optimisation | EPSI_Opt | 1
| Ultrafast 2D TOCSY | EPSI_Opt | 10–11
PSYCHE | Flip angle reference grid search | PSYCHE_FAScan | 301–331
| Optimisation (1 parameter) | PSYCHE_Opt | 101
| Optimisation (2 parameters) | PSYCHE_Opt | 102
| Optimisation (3 parameters) | PSYCHE_Opt | 103
| Optimisation (4 parameters) | PSYCHE_Opt | 104
| Pseudo-2D homodecoupled spectra | PSYCHE_Opt | 1400–1404
Solvent suppression | Offset reference grid search | Solvsupp_O1Scan | 101–141
| Optimisation (1 parameter) | Solvsupp_Opt | 11
| Optimisation (2 parameters) | Solvsupp_Opt | 14
| Optimisation (3 parameters) | Solvsupp_Opt | 15
| Optimisation (4 parameters) | Solvsupp_Opt | 16
| Representative optimised spectra | Solvsupp_Opt | 31, 35–37
DOSY | Optimisation (sequential) | DOSY_Opt | 101–102
| Optimised Oneshot DOSY | DOSY_Opt | 1, 1001–1021
| Optimisation (simultaneous) | DOSY_Opt | 201–202
Table S1: Datasets used for each section of the Supporting Information.

All of the POISE optimisation data was logged in poise.log files, which can be
found in the corresponding raw data folders. These log files detail the
optimisation trajectory and the cost function at each point, and can be parsed
using the parse_log() function (please see the documentation for more
information). However, some of these date from an older development version of
POISE in which less information was recorded, so some details (e.g. of the AU
programme used for acquisition) may be blank.

Full details of acquisition and processing parameters for each optimisation
can be found in these folders. For this reason, except for a few key values,
they are not tabulated in the main text or here.

###  S1.2 Figures

All figures are done in Python 3 and use the penguins package (written by
J.R.J.Y.) to parse and plot NMR data. As of the time of writing, this package
is still in development; consequently, the interface and functionality may
change in a backward-incompatible manner. In order to ensure full
reproducibility, please install version 0.4.1 of penguins, which was released
specifically for this paper. This can be done using the command

⬇

pip install penguins=0.4.1

After downloading the raw data, the NMR datasets themselves should be found in
a directory called datasets, and the scripts in a directory called figures.
The figures themselves are provided in this folder, but can be regenerated at
any point in time by re-running the scripts (the scripts also provide a
“trace” of which data is being used to plot the data).

The scripts will be executable as long as this directory structure is
preserved, as they use a relative path to locate the datasets. If the
directories are moved for any reason, the environment variable POISE_DATA
should be set such that it points to the datasets directory. For example, in
bash and similar shells, this can be done using (replace the path accordingly)

⬇

export POISE_DATA=/path/to/nmrpoise_data/datasets

##  S2 Routines

In POISE, routines are stored in the human-readable JSON format, and are thus
provided in the sections which follow. These JSON files should be placed in
the TopSpin directory .../exp/stan/nmr/
py/user/poise_backend/routines. Each routine consists of several entries:

  * •

name: The name of the routine, which can be arbitrarily chosen (although it
should match the filename: thus a routine with name my_routine should be
stored in my_routine.json).

  * •

pars: A list of string parameters which are to be optimised. These correspond
exactly to parameter names in TopSpin.

  * •

lb: A list of numbers acting as the lower bounds for each parameter. Units are
given in TopSpin’s native units, i.e. pulses are in microseconds, delays in
seconds, etc.

  * •

ub: A list of numbers acting as the upper bounds for each parameter.

  * •

init: A list of numbers acting as the initial values for each parameter.

  * •

tol: A list of numbers acting as the tolerances for each parameter.

  * •

cf: The name of the cost function. A function with this name must be defined
in either of the files poise_backend/costfunctions.py or
poise_backend/costfunctions_user.py. For a list of pre-installed cost
functions and descriptions, please refer to the POISE user documentation.

  * •

au: The name of the AU programme used for acquisition and processing. This can
be left blank (see below).

If no AU programme is specified by the user upon setup, the poise_1d or
poise_2d AU programmes will automatically be used depending on the
dimensionality of the dataset. These default AU programmes are used for many
of the optimisations described in the Supporting Information, and are detailed
below.

All the routines in the Supporting Information can be found in the
.../exp/stan/nmr/py/user/
poise_backend/example_routines directory, and can be directly copied to the
routines directory if desired.

###  S2.1 poise_1d AU programme

⬇

ZG // acquire

EFP // exponential apodisation, Fourier transform

APBK // automatic phase and baseline correction

QUIT

###  S2.2 poise_1d_noapk AU programme

This AU programme is not used by default, but comes installed together with
POISE, used in some of the optimisations described below.

⬇

ZG // acquire

EFP // exponential apodisation, Fourier transform

XCMD(”abs n”) // automatic baseline correction

QUIT

###  S2.3 poise_2d AU programme

⬇

ZG // acquire

XFB // 2D Fourier transform

XCMD(”apk2d”) // automatic phase correction

ABS2 // f2 baseline correction

QUIT

##  S3 1D ^1H spectra for all samples

As a reference, we provide here the ^1H spectra of all samples used in this
work. These can also be found in the raw data (Section S1).

Figure S1:  ^1H spectrum of $45\text{\,}\mathrm{m\textsc{m}}$ ferulic acid in
DMSO-$d_{6}$ ($500\text{\,}\mathrm{MHz}$ resonance frequency).  Figure S2:
^1H spectrum of $120\text{\,}\mathrm{m\textsc{m}}$ 3-fluorophenylboronic acid
in DMSO-$d_{6}$ ($700\text{\,}\mathrm{MHz}$ resonance frequency).  Figure S3:
^1H spectrum of $45\text{\,}\mathrm{m\textsc{m}}$ andrographolide in
DMSO-$d_{6}$ ($600\text{\,}\mathrm{MHz}$ resonance frequency).  Figure S4:
^1H spectrum of rodent urine in D2O ($400\text{\,}\mathrm{MHz}$ resonance
frequency). The water peak has been suppressed using the 1D NOESY /
presaturation sequence.

##  S4 $\mathbf{90^{\circ}}$ pulse width calibration

All spectra in this section were acquired on a $600\text{\,}\mathrm{MHz}$
Bruker AVIII spectrometer, equipped with a Prodigy N2 broadband cryoprobe with
a nominal maximum $z$-gradient strength of
$66\text{\,}\mathrm{G}\text{\,}{\mathrm{cm}}^{-1}$. The sample used was
$45\text{\,}\mathrm{m\textsc{m}}$ ferulic acid in DMSO-$d_{6}$.

###  S4.1 Cost function

The cost function used in this section, minabsint, comes pre-installed with
POISE. It can be mathematically expressed as

$f=\sum_{i}\sqrt{\mathbf{S}_{\text{real},i}^{2}+\mathbf{S}_{\text{imag},i}^{2}}$
(1)

where $\mathbf{S}_{\text{real}}$ and $\mathbf{S}_{\text{imag}}$ are vectors of
length SI corresponding to the real and imaginary parts of the spectrum
respectively. The elements $\mathbf{S}_{i}$ are the intensities of the
spectrum at each point along the spectral width.

###  S4.2 Reference grid search

The grid search, as used in the builtin popt, is a highly time-consuming
method of locating the optimum. It is used here only as a reference, i.e. to
prove that the faster optimisation algorithms are indeed finding the true
optimum.

Figure S5:  Dependence of the cost function on the pulse width (p1). The
minimum value of the cost function is $3.34\times 10^{7}$ and occurs at
$\mathtt{p1}=$48.3\text{\,}\mathrm{\SIUnitSymbolMicro s}$$.

###  S4.3 Optimisation

Entry | Method | Description | Optimum found ($\mathrm{\SIUnitSymbolMicro s}$) | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---|---
1 | popt |  MAGMIN cost function | 48.40 | 41 | 299
2 | pulsecal | – | 46.64 | – | 37
3a | POISE | Nelder–Mead | 48.38 | 10 | 76–79
3b |  | MDS | 48.38 | 10 | 77–80
3c |  | BOBYQA | 48.29–48.41 | 6–7 | 46–54
Table S2:  $None$ pulse width optimisations run with an initial point of
$48\text{\,}\mathrm{\SIUnitSymbolMicro s}$. popt optimisations were run
between values of $40\text{\,}\mathrm{\SIUnitSymbolMicro s}$ and
$56\text{\,}\mathrm{\SIUnitSymbolMicro s}$, with a linear increment of
$0.4\text{\,}\mathrm{\SIUnitSymbolMicro s}$. The POISE routine is: . The pulse
programme used was the Bruker standard zg. Key acquisition parameters: NS=1,
DS=0, SW=12.0, TD=16384, D1=1.5. $n_{\text{fev}}$ refers to the number of
function evaluations, i.e. the number of spectra acquired in the course of the
optimisation. The results are aggregated from five separate runs for each
optimisation algorithm.

The BOBYQA algorithm takes the exact value of the cost function into
consideration when deciding what points to sample. On the other hand, the
simplex-based methods only consider how the various points are ordered
according to the value of the cost function at each point. This is partly why
BOBYQA tends to converge faster than the other two algorithms, which in turn
explains why we have set it as the default algorithm for POISE; it also
accounts for the slight variability in the optimum found (note that this range
is still well within the specified tolerance of
$0.2\text{\,}\mathrm{\SIUnitSymbolMicro s}$, so it does not reflect a
failing). However, in edge cases when the cost function is extremely
unbalanced or noisy, we find that BOBYQA can be unduly biased by the resulting
cost function values (this is not the case here). We find that the simplex-
based NM and MDS methods, although slower on average, can be slightly more
reliable in such situations.

The next two tables describe optimisations starting from an inaccurate initial
point.

Entry | Method | Description | Optimum found ($\mathrm{\SIUnitSymbolMicro s}$) | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---|---
1a | POISE | Nelder–Mead | 48.38 | 14 | 109–114
1b |  | MDS | 48.38 | 14 | 108–112
1c |  | BOBYQA | 48.27–48.33 | 9 | 70
Table S3:  Pulse width optimisations with an initial point of $43\text{\,}\mathrm{\SIUnitSymbolMicro s}$. The POISE routine is the same as in Table S2, except with "init":[43.0]. All other acquisition parameters are identical. The results are aggregated from five separate runs for each optimisation algorithm.  Entry | Method | Description | Optimum found ($\mathrm{\SIUnitSymbolMicro s}$) | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---|---
1a | POISE | Nelder–Mead | 48.38 | 14 | 110–114
1b |  | MDS | 48.25–48.38 | 16 | 123–126
1c |  | BOBYQA | 48.26–48.33 | 9 | 69–70
Table S4:  Pulse width optimisations with an initial point of
$53\text{\,}\mathrm{\SIUnitSymbolMicro s}$. The POISE routine is the same as
in Table S2, except with "init":[53.0]. All other acquisittion parameters are
identical. The results are aggregated from five separate runs for each
optimisation algorithm.

As can be seen, POISE is capable of “recovering” from a bad start. Although
this requires more time, even in the worst-case scenario, the time taken is
still less than half of that needed for the full grid search (299 seconds, cf.
Entry 1, Table S2).

###  S4.4 Using POISE under automation

POISE provides a command-line interface which allows it to be called from
within AU and Python scripts once the appropriate routines and cost functions
have been set up. Thus, for example, a typical optimisation ‘procedure’ would
involve the following steps:

  1. 1.

Create a new dataset specially for the optimisation.

  2. 2.

Read in a parameter set, where the appropriate pulse programme for
optimisation, as well as all other unoptimised parameters, have been chosen.

  3. 3.

Run the POISE optimisation with poise <routine_name> <options>. (For the
available options please see the user documentation.)

  4. 4.

Transfer the optimised parameter values back to the original dataset, and run
the optimised experiment.

All of this can be easily automated using either AU or Python scripts. Here we
provide examples of both AU and Python scripts, called poisecal, which perform
a similar role to the existing pulsecal AU programme in TopSpin. Running
either of these will carry out a POISE optimisation in expno 99999, then set
the optimised value of p1 to the currently active dataset. (These examples are
also included in the documentation with extra explanation.) Please note that
both of these can be downloaded from GitHub; the links are near the top of the
source code.

In general, the AU programme should be preferred as it contains some extra
functionality (namely the option to dynamically specify the region of the
spectrum to optimise on: please see the header of the AU programme for more
information).

####  S4.4.1 AU programme

⬇

/******************************************************************************

* poisecal

* ——–

*

* Download from: https://git.io/JUHOk

*

* AU script which automates p1 optimisation using POISE. Behaves similarly to

* pulsecal: optimises the value of p1 in expno 99999, then sets the optimised

* value to the current dataset and shows a message.

*

* To optimise on only one particular part of the spectrum, either:

*

* 1) Use the command ‘dpl‘ to specify a particular region, then run

* ‘poisecal‘; or

* 2) use the command

*

* poisecal LEFT RIGHT

*

* where LEFT and RIGHT are the edges of the desired spectral window in

* ppm. The order in which these are specified does not matter, so

* ‘poisecal 1 4‘ and ‘poisecal 4 1‘ behave exactly the same way.

*

* To optimise only one end of a spectrum (e.g. above 4 ppm only), use option 2

* and pass a large value as the other argument, e.g. ‘poisecal 4 1000‘.

*

* SPDX-License-Identifier: GPL-3.0-or-later

*****************************************************************************/

GETCURDATA

int old_expno = expno;

double f1p, f2p;

// Argument parsing

if (i_argc == 4) { // argv[0] is always the programme name, argv[1] is ’exec’

double d1, d2;

char *cp;

d1 = strtod(i_argv[2], &cp);

if (*cp) STOPMSG(”Invalid arguments. Usage: poisecal [left right]”)

d2 = strtod(i_argv[3], &cp);

if (*cp) STOPMSG(”Invalid arguments. Usage: poisecal [left right]”)

f1p = d1 > d2 ? d1 : d2;

f2p = d1 > d2 ? d2 : d1;

}

else if (i_argc == 2) {

FETCHPAR(”F1P”, &f1p)

FETCHPAR(”F2P”, &f2p)

}

else STOPMSG(”Invalid arguments. Usage: poisecal [left right]”)

// Use EXPNO 99999 in the current folder for optimisation.

DATASET(name, 99999, procno, disk, user)

// Set some key parameters. Notice that these lines can be substantially cut

// if an appropriate parameter set is set up beforehand.

RPAR(”PROTON”, ”all”)

GETPROSOL

STOREPAR(”PULPROG”, ”zg”)

STOREPAR(”NS”, 1)

STOREPAR(”DS”, 0)

STOREPAR(”D 1”, 1.0)

STOREPAR(”RG”, 1)

STOREPAR(”F1P”, f1p)

STOREPAR(”F2P”, f2p)

STOREPARS(”F1P”, f1p)

STOREPARS(”F2P”, f2p)

// Run optimisation (uses BOBYQA by default).

XCMD(”sendgui xpy poise p1cal -q”)

// POISE stores the optimised value in p1 after it’s done. We can retrieve it

// here. Don’t try to get the *status* parameter, since that is not the

// optimised value (it is the value used for the last function evaluation!)

float p1opt;

FETCHPAR(”P 1”, &p1opt)

p1opt = p1opt/4;

// Move back to old dataset and set p1 to optimised value.

DATASET(name, old_expno, procno, disk, user)

VIEWDATA_SAMEWIN // not strictly necessary, just re-focuses the original
spectrum

STOREPAR(”P 1”, p1opt)

Proc_err(INFO_OPT, ”Optimised value of p1: %.3f”, p1opt);

// (Optional) Run acquisition.

// ZG

QUIT

// vim: ft=c

####  S4.4.2 Python script

⬇

”””

poisecalpy.py

————-

Download from: https://git.io/JUHOU

Python script which automates p1 optimisation using POISE. Behaves similarly
to

pulsecal: optimises the value of p1 in expno 99999, then sets the optimised

value to the current dataset and shows a message.

This does not have the command-line argument functionality of the poisecal AU

programme, and tends to suffer from some edge-case bugs which are not present

in poisecal. This is generally a result of TopSpin’s Python programmes being

more buggy than AU programmes. Nonetheless, as a result of this, we recommend

using the poisecal AU programme instead of this (which is retained mainly to

serve as a illustrative example in the documentation).

SPDX-License-Identifier: GPL-3.0-or-later

”””

# Read in F1P and F2P from current dataset.

f1p = GETPAR(”F1P”)

f2p = GETPAR(”F2P”)

# Use EXPNO 99999 in the current folder for optimisation.

old_dataset = CURDATA()

opt_dataset = CURDATA()

opt_dataset[1] = ”99999”

# Create new dataset and move to it.

NEWDATASET(opt_dataset, None, ”PROTON”)

RE(opt_dataset)

# Set some key parameters. Notice that these lines can be cut if an
appropriate

# parameter set is set up beforehand (and loaded using NEWDATASET()).

XCMD(”getprosol”)

PUTPAR(”PULPROG”, ”zg”)

PUTPAR(”NS”, ”1”)

PUTPAR(”DS”, ”0”)

PUTPAR(”D 1”, ”1”)

PUTPAR(”RG”, ”1”)

XCMD(”s f1p {}”.format(f1p)) # PUTPAR(”status F1P”) doesn’t work.

XCMD(”s f2p {}”.format(f2p)) # Even though the documentation says it should.

# Run optimisation (uses BOBYQA by default).

XCMD(”poise p1cal -q”)

# POISE stores the optimised value in p1 after it’s done. We can retrieve it

# here. Don’t try to get the *status* parameter, since that is not the

# optimised value (it is the value used for the last function evaluation!)

p1opt = float(GETPAR(”P 1”))/4

# Move back to old dataset and set p1 to optimised value.

RE(old_dataset)

PUTPAR(”P 1”, str(p1opt))

ERRMSG(”Optimised value of p1: {:.3f}”.format(p1opt))

# A TopSpin quirk: using MSG() will block subsequent commands until user hits

# ”OK”. You can use ERRMSG(), as is done here. There are other ways around

# this, see MSG_nonmodal() in the POISE frontend script.

# (Optional) Run acquisition.

# ZG()

##  S5 Ernst angle

All spectra in this section were acquired on a $500\text{\,}\mathrm{MHz}$
Bruker AVIII spectrometer, equipped with a triple-resonance TBO probe with a
nominal maximum $z$-gradient strength of
$50\text{\,}\mathrm{G}\text{\,}{\mathrm{cm}}^{-1}$. The sample used was
$45\text{\,}\mathrm{m\textsc{m}}$ ferulic acid in DMSO-$d_{6}$.

###  S5.1 Cost function

The cost function used in this section is maxrealint, which corresponds simply
to a summation of all peaks in the real part of the spectrum:

$f=-\sum_{i}\mathbf{S}_{\mathrm{real},i}$ (2)

The minus sign is necessary because optimisation algorithms always seek to
minimise the cost function. Thus, in this case, a larger peak intensity
corresponds to a more negative cost function, i.e. a “better” spectrum.

###  S5.2 Theoretical optima via $T_{1}$

In lieu of a reference grid search, we sought to measure the longitudinal
relaxation times ($T_{1}$) of each peak in the sample. The Ernst angle
$\theta$ which maximises the sensitivity for a given repetition time
$\tau_{\mathrm{r}}$ can then be calculated from this via

$\cos\theta=\exp\left(-\frac{\tau_{\mathrm{r}}}{T_{1}}\right)$ (3)

thus giving us the theoretical optimal value, without having to go through a
grid search. The $T_{1}$ values extracted from a gradient-enhanced inversion-
recovery experiment, and the corresponding Ernst angles (for a repetition time
of $\tau_{\mathrm{r}}=$1.20\text{\,}\mathrm{s}$$, as used in the
optimisation), are summarised in Table S5.

Peak |  ^1H chemical shift (ppm) |  $T_{1}$ ($\mathrm{s}$) |  $\theta$ (∘)
---|---|---|---
1 | 7.49 | 1.750 | 59.8
2 | 7.27 | 0.977 | 73.0
3 | 7.08 | 1.279 | 67.0
4 | 6.79 | 1.615 | 61.6
5 | 6.36 | 1.415 | 64.6
6 | 3.49 | 0.949 | 73.6
Table S5:  $T_{1}$ (at $500\text{\,}\mathrm{MHz}$) and Ernst angles for each
peak in ferulic acid, calculated for a repetition time of
$1.20\text{\,}\mathrm{s}$.

###  S5.3 Optimisation

It is possible to show, both theoretically and experimentally, that in order
to maximise the sensitivity per unit time, the relaxation delay should be set
to zero and the Ernst angle applied with $\tau_{\mathrm{r}}$ simply equal to
the acquisition time AQ. See, e.g., Concepts Magn. Reson. 1992, 4 (2), 153–160
(DOI: 10.1002/cmr.1820040204) and J. Mol. Spectrosc. 1970, 35 (2), 298–305
(DOI: 10.1016/0022-2852(70)90205-5).

Furthermore, in order to make an Ernst angle optimisation worthwhile, the
repetition time must be on the order of $T_{1}$. With longer repetition times
on the order of $3T_{1}$ or so, the Ernst angle is very close to $None$, which
transforms the issue into the $None$ calibration already described in the
previous section. In the event, we elected to use a repetition time of
$1.20\text{\,}\mathrm{s}$, which is on the order of the $T_{1}$ times shown in
Table S5.

In the first optimisation, we set the region under optimisation to be just the
aromatic and olefinic peaks between 6 and $8\text{\,}\mathrm{ppm}$ (i.e. peaks
1–5). Instead of changing the actual spectral window, this is most easily done
using the dpl command in TopSpin; for more information, please refer to the
user documentation.

Entry | Algorithm | Optimum found (∘) | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---
1 | Nelder–Mead | 67.5–73.1 | 9–13 | 91–132
2 | MDS | 67.5–73.1 | 9 | 90–92
3 | BOBYQA | 70.1–70.7 | 7 | 70–71
Table S6:  Ernst angle optimisations on all aromatic and olefinic peaks in
ferulic acid (between 6 and $8\text{\,}\mathrm{ppm}$). The POISE routine used
here is: . The pulse programme used was zgvfa_jy (Section S5.4). Key
acquisition parameters: NS=1, DS=4, D1=0, SW=14, TD=16384, AQ=1.17.
$n_{\text{fev}}$ refers to the number of function evaluations, i.e. the number
of spectra acquired in the course of the optimisation. The results are
aggregated from five separate runs for each optimisation algorithm.

Of course, the user may select the spectral bounds themselves in order to
maximise the overall sensitivity for the entire spectrum, a particular region,
or even a specific peak of interest. In the second case, we ran an
optimisation on only the peak at $6.79\text{\,}\mathrm{ppm}$.

Entry | Algorithm | Optimum found ($\mathrm{s}$) | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---
1 | Nelder–Mead | 60.0–67.5 | 9–11 | 91–111
2 | MDS | 65.6–67.5 | 11 | 110–111
3 | BOBYQA | 60.0–65.2 | 6–7 | 59–71
Table S7:  Ernst angle optimisations on the peak at
$6.79\text{\,}\mathrm{ppm}$ in ferulic acid. The POISE routine and all
acquisition parameters are the same as in Table S6, but the spectral region
under optimisation was set to be 6.71–$6.87\text{\,}\mathrm{ppm}$ via the dpl
TopSpin command. $n_{\text{fev}}$ refers to the number of function
evaluations, i.e. the number of spectra acquired in the course of the
optimisation. The results are aggregated from five separate runs for each
optimisation algorithm.

###  S5.4 Variable flip angle pulse-acquire pulse programme

Note that the following pulse programme, which requires p1 to be the
calibrated $None$ pulse width and sets cnst20 to the flip angle, should only
be used if the actual flip angle is of interest.

If the sole aim is to maximise sensitivity, then the standard Bruker zg pulse
programme suffices. The routine should be modified to optimise p1 directly
instead of the flip angle itself. In essence, this means that we do not care
what the flip angle is, we merely want the pulse width which maximises the
sensitivity. This approach has an advantage in that p1 does not need to be
calibrated beforehand.

⬇

;zgvfa_jy

;

;download from: https://git.io/JnokS

;

;pulse-acquire with variable flip angle (set by cnst20)

;set p1 to be calibrated 90 degree pulse width

;Jonathan Yong - University of Oxford - 20 June 2021

;

;$CLASS=HighRes

;$DIM=1D

;$TYPE=

;$SUBTYPE=

;$COMMENT=

#include <Avance.incl>

”p20=p1*cnst20/90”

”acqt0=-p20*2/3.1416”

1 ze

2 30m

d1

p20 ph1

go=2 ph31

30m mc #0 to 2 F0(zd)

exit

ph1=0 2 2 0 1 3 3 1

ph31=0 2 2 0 1 3 3 1

;pl1 : f1 channel - power level for pulse (default)

;p1 : f1 channel - high power pulse

;d1 : relaxation delay; 1-5 * T1

;ns: 1 * n, total number of scans: NS * TD0

;cnst20: excitation flip angle in degrees

##  S6 NOE mixing time

All spectra in this section were acquired on a $700\text{\,}\mathrm{MHz}$
Bruker AVIII spectrometer, equipped with a TCI H/C/N cryoprobe with a nominal
maximum $z$-gradient strength of
$53\text{\,}\mathrm{G}\text{\,}{\mathrm{cm}}^{-1}$. The sample used was
$120\text{\,}\mathrm{m\textsc{m}}$ 3-fluorophenylboronic acid in DMSO-$d_{6}$.

###  S6.1 Cost function

The cost function used for the optimisation here, named noe_1d, is pre-
installed with POISE. It is only applicable to 1D selective NOE spectra. In
short, it takes the processed spectrum, removes a region of
$50\text{\,}\mathrm{Hz}$ around the selectively excited peak, and calculates
the intensity of the remaining parts of the spectrum. Since we aim to maximise
this intensity, the cost function must return the negative of this intensity.
Further details are best explained through the comments in the source code,
which is reproduced here:

⬇

”””

costfunctions.py

—————-

Contains default cost functions. Custom cost functions should be added in

costfunctions_user.py (using the same format).

You can add custom cost functions in here (they will work), but they will be

shadowed by any similarly named cost functions in costfunctions_user.py. Also,

this file will be overwritten if POISE is reinstalled.

SPDX-License-Identifier: GPL-3.0-or-later

”””

import numpy as np

from .cfhelpers import *

from .shared import _g

def noe_1d():

”””

Measures the intensity of peaks in the spectrum, *except* for anything

within 25 Hz of the selectively excited peak.

”””

# SPOFFS2 plus O1 is the frequency of the selective pulse, in Hz.

# We can get the values of both using getpar(), which automatically

# converts numeric parameters to floats.

f = getpar(”SPOFFS2”) + getpar(”O1”)

# Define the (rough) bandwidth of the selective pulse. We will ignore the

# region of the spectrum that is bw Hz wide and centred on the selective

# pulse.

bw = 50.0

# To get the spectrum *without* the ignored region, we split up the

# spectrum into two ”portions”. We need to construct the ’bounds’ parameter

# for both portions.

sfo1 = getpar(”SFO1”)

lowerhalf = f”{(f + bw/2)/sfo1:.3f}..” # frequencies (f + bw)/2 and above

upperhalf = f”..{(f - bw/2)/sfo1:.3f}” # frequencies up to (f - bw)/2

# We use get1d_real() to get both portions of the spectrum, passing the

# bounds as appropriate. The intensities of the spectrum are obtained by

# summing all points (np.sum).

# Note that manually passing the bounds parameter will override the F1P/F2P

# parameters if set by the user (this is probably a good thing).

upper_integral = np.sum(get1d_real(bounds=upperhalf))

lower_integral = np.sum(get1d_real(bounds=lowerhalf))

# Lastly, we take the negative absolute value of the resulting sum, in

# order to make sure that the resulting cost function is negative

# regardless of how the NOE peaks are phased (positive or negative).

return -abs(upper_integral + lower_integral)

def minabsint():

”””

Cost function which minimises the absolute (magnitude-mode) intensity of

the spectrum. This is probably the easiest cost function. :-)

”””

# Get the real part of the spectrum as a numpy ndarray. Because no bounds

# are explicitly passed, this defaults to the whole spectrum if F1P/F2P are

# not specified. If F1P/F2P *are* specified (either via the ‘dpl‘ TopSpin

# command, or manually), this will only return the part of the spectrum

# between those two values.

r = get1d_real()

# …and the imaginary part

i = get1d_imag()

# The magnitude-mode spectrum is obtained by combining these.

mag = np.abs(r + 1j * i)

# The intensity of the magnitude-mode spectrum is just the sum of all

# points.

return np.sum(mag)

def maxabsint():

”””

Cost function which maximises the absolute (magnitude-mode) intensity of

the spectrum.

”””

r = get1d_real()

i = get1d_imag()

# This is the same as minabsint except that we have a negative sign.

# Because the optimisation always seeks to *minimise* the cost function,

# this essentially tries to *maximise* np.sum(…), i.e. maximise the

# spectral intensity.

return -np.sum(np.abs(r + 1j*i))

def minrealint():

”””

Minimises the intensity of the real part of the spectrum. If the spectrum

has negative peaks this cost function will try to maximise those.

”””

return np.sum(get1d_real())

def maxrealint():

”””

Maximises the intensity of the real part of the spectrum.

”””

return -np.sum(get1d_real())

def zeronetrealint():

”””

Tries to get the *net* intensity of the real spectrum to be as close to

zero as possible. Note that dispersion-mode peaks will contribute to this

cost function (as they add to zero).

”””

return np.abs(np.sum(get1d_real()))

def zerorealint():

”””

Tries to get the real spectrum to be as close to zero as possible. This

takes the absolute value of the real spectrum before summing, so

dispersion-mode peaks do not add to zero.

”””

return np.sum(np.abs(get1d_real()))

def zerorealint_squared():

”””

Tries to get the real spectrum to be as close to zero as possible. This

uses a sum-of-squares metric, which (compared to “zerorealint“) leads to

two main differences. Firstly, dispersion-mode peaks do not sum to zero

because sign information is lost upon squaring. Secondly, it more strongly

penalises large deviations away from zero.

”””

s = get1d_real()

return np.sum(s * s)

def epsi_gradient_drift():

”””

Calculates the amount of ’gradient drift’ seen in a 1D EPSI acquisition, as

reflected by the position of the ’echo’ moving over time. This can be

caused by imbalanced positive and negative gradients.

”””

from numpy.polynomial.polynomial import Polynomial

# — Calculate key parameters ——————————-

fid = get1d_fid(remove_grpdly=True)

# Number of complex points in k-space, i.e. number of points per EPSI

# gradient. TD2 is total number of real & imag points in FID. L3 is number

# of EPSI loops (one loop includes both pos + neg gradient). So TD2 / (2 *

# L3) is the number of points in one EPSI gradient (i.e. only positive

# gradient). The extra factor of 2 is because we’re interested only in

# complex points.

td_k = int(getpar(”TD”) / (2 * 2 * getpar(”L3”)))

# Time between consecutive EPSI positive gradients, in seconds

td_t2 = int(getpar(”L3”))

dw_eff = getpar(”AQ”) / td_t2

# — Perform EPSI processing ——————————–

# Reshape into 2D matrix

ser = fid.reshape((-1, td_k))

# Discard the part acquired with negative gradients

ser = ser[0::2, :]

# Discard any part of the spectrum that was not acquired during an EPSI

# gradient, i.e. if the delay D6 was nonzero.

if getpar(”D6”) > 0:

ineligible_epsi_points = int(1e6 * getpar(”D6”) / (getpar(”DW”) * 2))

td_k = td_k - ineligible_epsi_points

ser = ser[:, :td_k]

# — Apodisation ——————————————–

# Along k-dimension (Hamming window)

alpha_0 = 0.54

ks = np.linspace(0, 1, td_k)

k_winfunc = (alpha_0

- (1 - alpha_0) * np.cos(2 * np.pi * np.linspace(0, 1, td_k)))

ser = ser * k_winfunc[np.newaxis, :]

# Along direct dimension

t2_winfunc = np.sin(np.pi * np.linspace(0, 1, td_t2))

ser = ser * t2_winfunc[:, np.newaxis]

abs_ser = np.abs(ser)

# Calculate k- and t2-axes

t2_values = np.arange(td_t2) * dw_eff

k_values = np.linspace(-0.5, 0.5, td_k)

# Drop all rows (i.e. all values of t2) for which the maximum is less than

# 20% of the overall maximum.

threshold_amp = 0.2 * np.max(abs_ser)

maxima_along_rows = np.max(abs_ser, axis=1)

indices_to_use = np.nonzero(maxima_along_rows > threshold_amp)

# Locate the maxima

maximal_indices_along_k = np.argmax(abs_ser, axis=1)

# Get the value of k at each maximum

k_max_values = k_values[maximal_indices_along_k]

poly = Polynomial.fit(x=t2_values[indices_to_use],

y=k_max_values[indices_to_use], deg=1)

return np.abs(poly.coef[1])

# The cost functions below were used in the POISE paper, but are less likely
to

# be generally useful. If you want to use them, simply uncomment them (remove

# the three single quotes around them).

”’

def asaphsqc():

”””

Maximises the intensity of the projection of a 2D spectrum onto the f2

axis.

”””

# First we get the ’rr’ part of the spectrum (the displayed part). This

# assumes that it has been phased already (by the AU program).

# This is a 2D numpy array:

# [ [1 , 2, 3, …, si2],

# [. , ., ., …, …],

# [. , ., ., …, …],

# [si1, ., ., …, …]]

# where si1 and si2 are the values of SI in the indirect and direct

# dimension respectively.

spec = get2d_rr()

# We take the ’skyline’ projection, i.e. the highest point along each

# column. This is basically equivalent to TopSpin’s ‘f2projp‘ command.

# If you want a projection onto f1 instead of f2, you can use axis=1.

proj = np.amax(spec, axis=0)

return -np.sum(proj)

”’

”’

def specdiff():

”””

Returns the ”spectral difference” between the current spectrum and a

reference spectrum with EXPNO 1 and PROCNO 1\.

”””

# We can construct the path to the reference spectrum using make_p_spec().

reference_path = make_p_spec(expno=1, procno=1)

# The way this is currently coded, the reference spectrum has to be placed

# in EXPNO 1 and PROCNO 1. In theory, this can be specified using a TopSpin

# parameter in the optimisation dataset, for example ”cnst30” (or any other

# random cnst). In order to use this, just add:

# expno = # int(getpar(”cnst30”))

# and then after that you can use

# make_p_spec(expno=expno, procno=1)

# Then we can get the reference, or ’target’, spectrum as a numpy array.

# Note that this will respect the F1P/F2P parameters, *not* of the

# reference spectrum, but of the spectrum being optimised. This is because

# get1d_real() does not actually read F1P and F2P from the spectrum. It

# just uses the global values stored in _g.spec_f1p and _g.spec_f2p, which

# refer to the F1P/F2P values of the spectrum being optimised.

target = get1d_real(p_spec=reference_path)

# … and the current spectrum that is being optimised.

spec = get1d_real()

# Finally, return the 2-norm of the difference between the two normalised

# spectra.

return np.linalg.norm(target/np.linalg.norm(target) \-

spec/np.linalg.norm(spec))

”’

”’

def dosy():

”””

Tries to make a spectrum (1D diffusion experiment with variable gradient

strength) have 25% intensity of a target spectrum (1D diffusion experiment

with minimum gradient strength). This should be used in conjunction with

the dosy_opt Python programme (after POISE is installed via ‘pip‘, this can

be installed from the TopSpin command-line via ‘poise –install dosy‘.)

”””

# The target spectrum is hardcoded to have EXPNO 99998. See psyche() for

# more discussion about this.

# As in psyche(), the F1P/F2P parameters are not taken from the reference

# spectrum but from the spectrum being optimised.

reference_path = make_p_spec(expno=99998, procno=1)

target = get1d_real(p_spec=reference_path)

# This is the spectrum being optimised.

spec = get1d_real()

# The intensity of the optimised spectrum is np.sum(spec), and likewise for

# the target spectrum. We take the ratio of the two and measure how ”far

# away” this is from the ideal value of 0.25 (which corresponds to 75%

# attenuation).

return np.abs(np.sum(spec)/np.sum(target) - 0.25)

”’

”’

def dosy_aux():

”””

Non-absolute value of dosy(). To be used in the first stage of dosy_opt.

”””

reference_path = make_p_spec(expno=99998, procno=1)

target = get1d_real(p_spec=reference_path)

spec = get1d_real()

return np.sum(spec)/np.sum(target) - 0.25

”’

”’

def dosy_2p():

”””

Two-parameter DOSY optimisation.

”””

reference_path = make_p_spec(expno=99998, procno=1)

target = get1d_real(p_spec=reference_path)

spec = get1d_real()

return np.abs(np.sum(spec)/np.sum(target) - 0.25) + getpar(”D20”)

”’

###  S6.2 Reference grid search

The grid search was carried out using actual 2D NOESY spectra; this has the
advantage of proving that the mixing time optimised using 1D spectra does in
fact accurately represent the best mixing time for 2D spectra. The intensity
of each crosspeak here refers to the intensity of the most negative point
within a $0.12\text{\,}\mathrm{ppm}$ by $0.12\text{\,}\mathrm{ppm}$ box
centred on the crosspeak. A broad minimum around a mixing time of 3–4 seconds
can be seen.

Figure S6: Dependence of crosspeak intensity on the NOE mixing time (d8).

###  S6.3 Optimisation

Entry | Algorithm | Optimum found ($\mathrm{s}$) | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---
1 | Nelder–Mead | 3.25–3.88 | 16–18 | 268–312
2 | MDS | 3.63–3.75 | 16–18 | 269–305
3 | BOBYQA | 3.38–3.80 | 6–10 | 88–162
Table S8:  NOE mixing time optimisations. The POISE routine used here is: .
The pulse programme used here is spfgsenoezs2.1_jy (Section S6.4). Key
acquisition parameters: NS=2, DS=1, D1=1. $n_{\text{fev}}$ refers to the
number of function evaluations, i.e. the number of spectra acquired in the
course of the optimisation. The results are aggregated from five separate runs
for each optimisation algorithm.

###  S6.4 1D NOE pulse programme

The 1D NOE spectra were acquired with a modified pulse programme with two
inversion pulses during the NOE mixing time to minimise artefacts arising from
relaxation. (In principle, these artefacts do not contribute to the cost
function as they ideally sum to zero due to their antiphase nature. However,
it is almost always a good idea to try to minimise noise in the cost function
where possible.) The pulse programme is attached below, and can also be found
in the main POISE repository. For most compounds, the existing family of
Bruker pulse sequences selnogpzs should suffice.

⬇

;spfgsenoezs2.1_jy

;

;download from: https://git.io/JUE0P

;

;Tim Claridge & Jonathan Yong - University of Oxford

;1D NOESY using selective refocussing with two shaped pulses

;dipolar coupling may be due to noe or chemical exchange.

;H. Kessler, H. Oschkinat, C. Griesinger & W. Bermel,

; J. Magn. Reson. 70, 106 (1986)

;J. Stonehouse, P. Adell, J. Keeler & A.J. Shaka, J. Am. Chem. Soc 116,

; 6037 (1994)

;K. Stott, J. Stonehouse, J. Keeler, T.L. Hwang & A.J. Shaka,

; J. Am. Chem. Soc 117, 4199-4200 (1995)

;Tim C 2/2008 Using two nulling 180s and ZQC filter at end of mixing time

;Includes adiabatic pulse/gradient combination for ZQC suppression

;Thrippleton & Keeler, Angew. Chem. Int. Ed., 2003, 42, 3938-3941.

;Tim C 10/2007 updated to match Bruker format

;Tim C 11/2008 matches Bruker ZQS nomenclature

;Jon Y 21/07/2020 added CALC_SPOFFS acquisition flag

;$CLASS=HighRes

;$DIM=1D

;$TYPE=

;$SUBTYPE=

;$COMMENT=

#include <Avance.incl>

#include <Grad.incl>

”p2=p1*2”

”d21=d8*0.3-p16-d16-3u”

”d22=d8*0.5-p16*2-d16*2-6u”

”d23=d8*0.2-p16-d16*2-p32-53u”

#ifdef CALC_SPOFFS

”spoff2=bf1*(cnst21/1000000)-o1”

#else

#endif /*CALC_SPOFFS*/

1 ze

2 30m

20u BLKGRAD

d1

50u UNBLKGRAD

(p1 ph1):f1

3u

p16:gp1

d16 pl0:f1

3u

p12:sp2:f1 ph2:r

3u

p16:gp1

d16 pl1:f1

3u

(p1 ph4):f1

d21

p16:gp3

d16

3u

(p2 ph5):f1

3u

p16:gp3*-1

d16

d22

p16:gp4

d16

3u

(p2 ph5):f1

3u

p16:gp4*-1

d16

d23

20u pl0:f1

10u gron0 ;ZQC suppression

(p32:sp29 ph7):f1

20u groff

d16 pl1:f1

(p1 ph6):f1

go=2 ph31

30m mc #0 to 2 F0(zd)

20u BLKGRAD

exit

ph1=0 2

ph2=0 0 1 1 2 2 3 3

ph4=0

ph5=0

ph6=0

ph7=0

ph31=0 2 2 0 0 2 2 0

;pl0 : 120dB

;pl1 : f1 channel - power level for pulse (default)

;sp2: f1 channel - shaped pulse

;sp29: f1 channel - shaped pulse (adiabatic sweep)

;p1 : f1 channel - 90 degree high power pulse

;p2 : f1 channel - 180 degree high power pulse

;p12: f1 channel - 180 degree shaped pulse

;p16: homospoil/gradient pulse [1 msec]

;p32: f1 channel - 180 shaped adiabatic pulse

; smoothed CHIRP

;d1 : relaxation delay; 1-5 * T1

;d8 : mixing time

;d16: delay for homospoil/gradient recovery

;d21: d8*0.3-p16-d16-3u

;d22: d8*0.5-p16*2-d16*2-6u

;d23: d8*0.2-p16-d16*2-p32-53u

;NS: 2 * n

;DS: 4

;cnst21: chemical shift for selective pulse (offset, in ppm)

;choose p12 according to desired selectivity

;the flip-angle is determined by the amplitude

;set O1 on resonance on the multiplet to be excited or use spoffs

;use gradient ratio: gp 0 : gp 1 : gp 3 : gp 4

; ~10 : 17 : 40 : 25

;for z-only gradients:

;gpz0: 10% or so

;gpz1: 17%

;gpz3: 40%

;gpz4: 25%

;use gradient files:

;gpnam1: SINE.100

;gpnam2: SINE.100

;gpnam3: SINE.100

;gpnam4: SINE.100

##  S7 Inversion-recovery

All spectra in this section were acquired on a $500\text{\,}\mathrm{MHz}$
Bruker AVIII spectrometer, equipped with a triple-resonance TBO probe with a
nominal maximum $z$-gradient strength of
$50\text{\,}\mathrm{G}\text{\,}{\mathrm{cm}}^{-1}$. The sample used was
$45\text{\,}\mathrm{m\textsc{m}}$ ferulic acid in DMSO-$d_{6}$.

###  S7.1 Cost function

In this optimisation, we sought to make the real part of the spectrum as close
to zero as possible: this indicates the “null” in the inversion-recovery
experiment where the delay $\tau$ is exactly $T_{1}\ln 2$. This is
accomplished by the zerorealint cost function, which comes installed with
POISE:

$f=\sum_{i}\left\lvert\mathbf{S}_{\mathrm{real},i}\right\rvert$ (4)

###  S7.2 2D inversion-recovery

$T_{1}$ for each peak was measured via a gradient-selected inversion-recovery
experiment. The data presented here is the same as in Table S5, but
additionally includes the theoretical optimum of $T_{1}\ln 2$ for each peak.

Peak |  ^1H chemical shift (ppm) |  $T_{1}$ ($\mathrm{s}$) |  $T_{1}\ln 2$ ($\mathrm{s}$)
---|---|---|---
1 | 7.49 | 1.750 | 1.213
2 | 7.27 | 0.977 | 0.677
3 | 7.08 | 1.279 | 0.887
4 | 6.79 | 1.615 | 1.119
5 | 6.36 | 1.415 | 0.981
6 | 3.49 | 0.949 | 0.658
Table S9:  $T_{1}$ (at $500\text{\,}\mathrm{MHz}$) and the corresponding
theoretical null in the inversion-recovery profile, given by $T_{1}\ln 2$, for
each peak in ferulic acid.

###  S7.3 Optimisation

As previously mentioned, it is possible to use the dpl TopSpin command to
restrict the spectral region that the cost function acts on. In the first
optimisation, we set this region to be between 6 and $8\text{\,}\mathrm{ppm}$.
This provides an “averaged” value of $T_{1}\ln 2$, which can be useful in e.g.
choosing relaxation delays for 2D experiments, which invariably need to be a
compromise between multiple peaks.

Entry | Algorithm | Optimum found ($\mathrm{s}$) | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---
1 | Nelder–Mead | 0.938–0.969 | 14–16 | 204–235
2 | MDS | 0.956–0.975 | 16 | 233–235
3 | BOBYQA | 0.953–0.971 | 9–11 | 130–160
Table S10:  Inversion-recovery optimisations on all aromatic and olefinic
peaks in ferulic acid (between 6 and $8\text{\,}\mathrm{ppm}$. The POISE
routine used here is: . The pulse programme used was t1irgp1d_jy (Section
S7.4). Key acquisition parameters: NS=1, DS=0, D1=5, AQ=3.6 (note that there
is also an additional ca. 5 seconds between experiments due to pulse programme
compilation, etc.). $n_{\text{fev}}$ refers to the number of function
evaluations, i.e. the number of spectra acquired in the course of the
optimisation. The results are aggregated from five separate runs for each
optimisation algorithm.

In the second optimisation we focused on only the peak at
$7.08\text{\,}\mathrm{ppm}$. This provides a relatively “fast” method of
determining $T_{1}$ for one specific peak in a spectrum.

Entry | Algorithm | Optimum found ($\mathrm{s}$) | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---
1 | Nelder–Mead | 0.863–0.875 | 14 | 202–205
2 | MDS | 0.863–0.869 | 14 | 203–204
3 | BOBYQA | 0.862–0.873 | 9–10 | 128–145
Table S11:  Inversion-recovery optimisations on only the peak at
$7.08\text{\,}\mathrm{ppm}$. The POISE routine and all acquisition parameters
are the same as in Table S10, but the spectral region under optimisation was
set to be 7.02–$7.15\text{\,}\mathrm{ppm}$ via the dpl TopSpin command.
$n_{\text{fev}}$ refers to the number of function evaluations, i.e. the number
of spectra acquired in the course of the optimisation. The results are
aggregated from five separate runs for each optimisation algorithm.

###  S7.4 1D inversion-recovery pulse programme

⬇

;t1irgp1d_jy

;

;download from: https://git.io/Jnok6

;

;T1 measurement using inversion recovery

;using broadband 1H inversion pulse and field gradients

;1D version with only one increment

;set inversion-recovery delay to d27

;

;Tim Claridge & Jonathan Yong - University of Oxford - 20 June 2021

;

;$CLASS=HighRes

;$DIM=1D

;$TYPE=

;$SUBTYPE=

;$COMMENT=

#include <Avance.incl>

#include <Delay.incl>

#include <Grad.incl>

”acqt0=-p1*2/3.1416”

1 ze

2 50m BLKGRAD

d1

”TAU=(d27-p19-d16-4u)”

if ”TAU > 0.96”

{

”DELTA=TAU-0.45”

}

50u UNBLKGRAD

4u pl0:f1

(p44:sp30 ph1):f1

p19:gp1

d16

if ”TAU > 0.96”

{

0.4s

50m BLKGRAD

DELTA pl1:f1

}

else

{

TAU pl1:f1

}

p1 ph2

go=2 ph31

50m mc #0 to 2 F0(zd)

50m BLKGRAD

exit

ph1=0 2

ph2=0 0 2 2 1 1 3 3

ph31=0 0 2 2 1 1 3 3

;pl1 : f1 channel - power level for pulse (default)

;p1 : f1 channel - 90 degree high power pulse

;d1 : relaxation delay; 1-5 * T1

;ns: 1

;ds: 0

;spnam30: 1H inversion [Bip720,50,20.1]

;p44: broadband inversion pulse [240 us]

;d27: inversion-recovery delay

##  S8 ASAP-HSQC

All spectra in this section were acquired on a $700\text{\,}\mathrm{MHz}$
Bruker AVIII spectrometer, equipped with a TCI H/C/N cryoprobe with a nominal
maximum $z$-gradient strength of
$53\text{\,}\mathrm{G}\text{\,}{\mathrm{cm}}^{-1}$. The sample used was
$120\text{\,}\mathrm{m\textsc{m}}$ 3-fluorophenylboronic acid in DMSO-$d_{6}$.

###  S8.1 Cost function

The cost function used here, asaphsqc, simply takes the $f_{2}$ projection of
the acquired 2D spectrum and calculates the total intensity by summing up all
points. Unlike the previous cost functions, it is not builtin: to enable this
cost function, please uncomment the asaphsqc() function in
poise_backend/costfunctions.py.

###  S8.2 Reference grid search

Crosspeak intensities here refer to the maximum peak height, i.e. the
intensity of the highest point found within a $3\text{\,}\mathrm{ppm}$
($f_{1}$) by $0.12\text{\,}\mathrm{ppm}$ ($f_{2}$) box centred on the
crosspeak.

Figure S7: Dependence of crosspeak intensity on the INEPT delay
$\Delta=(1/4\cdot\mathtt{cnst3})$.

The exact point of maximum intensity is not obvious (and not really useful in
this case, as essentially any point near the maximum will give essentially the
same result). However, broadly speaking, it can be seen that the intensity
reaches a maximum at around $230\text{\,}\mathrm{H}\mathrm{z}$, after which it
plateaus off.

###  S8.3 Optimisation

Entry | Algorithm | Optimum found | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---
1 | Nelder–Mead | 250.0–256.3 | 8–9 | 169–195
2 | MDS | 237.5–243.8 | 8 | 171–179
3 | BOBYQA | 229.8–245.6 | 4–7 | 114–157
Table S12:  ASAP-HSQC INEPT delay optimisations. The POISE routine used was: .
The pulse programme used here is asaphsqc_jy (Section S8.4). Key acquisition
parameters: NS=1, DS=2, TD1=36, TD2=512, D1=0.1. $n_{\text{fev}}$ refers to
the number of function evaluations, i.e. the number of spectra acquired in the
course of the optimisation. The results are aggregated from five separate runs
for each optimisation algorithm.

###  S8.4 ASAP-HSQC pulse programme

The pulse programme used here is almost identical to the published version but
includes several minor changes for convenience.

⬇

;asaphsqc_jy

;

;download from: https://git.io/JUHOo

;

;modified from version by Luy et al. to include internal power calculations
for

;DIPSI-2 mixing, as well as URPs

;order of gradients also switched round to correct frequencies

;added extra delay d30 before experiment starts

;Jonathan Yong - University of Oxford - 26 August 2020

;

;ASAP-HSQC

;2D H-1/X correlation via double inept transfer

;phase sensitive using Echo/Antiecho-TPPI gradient selection

;with decoupling during acquisition

;with DIPSI-2 mixing after acquisition for ASAP

;using shaped pulses on f1 and f2.

;David Schulze-Sunninghausen, Johanna Becker, Burkhard Luy

;”Rapid Heteronuclear Single Quantum Correlation NMR Spectra at Natural
Abundance”

;J. Am. Chem. Soc. 2014, 136 (4), 1242-1245

#include <Avance.incl>

#include <Grad.incl>

#include <Delay.incl>

”cnst12=1/(1+sfo2/sfo1)”

”cnst13=(1-sfo2/sfo1)”

define list<gradient> EA1 = { 1 cnst12 }

define list<gradient> EA2 = { cnst13 1 }

”p2=p1*2”

”d0=3u”

”d11=30m”

”d16=200u”

”d20=1m”

”FACTOR1=(d9/(p6*115.112))/2”

”l1=FACTOR1*2” ; TOCSY loop counter

”d0=3u”

”in0=inf1/2”

”d4=1s/(4*cnst4)” ; 1/(4J)XH

”d3=1s/(4*cnst3)” ; reduced d4

”DELTA=0.5*(2*d3+p21-p1)”

”DELTA1=0.5*(2*d3-p21+p1)”

”DELTA4=d4-p16-d16-3u”

”DELTA3=p16+3u+d16+p21-p25+p1”

”DELTA2=0.5*(DELTA3+p16+d16+3u-2*d0)”

”cnst21=20000”

”spw21=plw1*(cnst21/(1000000/(4*p1)))*(cnst21/(1000000/(4*p1)))”

”spw22=spw21”

”cnst23=10000”

”spw23=plw2*(cnst23/(1000000/(4*p3)))*(cnst23/(1000000/(4*p3)))”

”spw24=spw23”

”spw25=spw23”

”spw26=spw23”

; calculate mixing length

”cnst9=p6*115.112*l1/1000”

”cnst9=cnst9” ; so that it’s displayed in ased

; calculate 90deg pulse duration & power level for mixing

”p6=(1000000/(4*cnst10))”

”plw10=plw1*(p1/p6)*(p1/p6)”

1 ze

d30

d11 pl12:f2

2 d1 do:f2

3 d20 UNBLKGRAD

p16:gp2

d16 pl10:f1

;begin DIPSI2

4 p6*3.556 ph23

p6*4.556 ph25

p6*3.222 ph23

p6*3.167 ph25

p6*0.333 ph23

p6*2.722 ph25

p6*4.167 ph23

p6*2.944 ph25

p6*4.111 ph23

p6*3.556 ph25

p6*4.556 ph23

p6*3.222 ph25

p6*3.167 ph23

p6*0.333 ph25

p6*2.722 ph23

p6*4.167 ph25

p6*2.944 ph23

p6*4.111 ph25

p6*3.556 ph25

p6*4.556 ph23

p6*3.222 ph25

p6*3.167 ph23

p6*0.333 ph25

p6*2.722 ph23

p6*4.167 ph25

p6*2.944 ph23

p6*4.111 ph25

p6*3.556 ph23

p6*4.556 ph25

p6*3.222 ph23

p6*3.167 ph25

p6*0.333 ph23

p6*2.722 ph25

p6*4.167 ph23

p6*2.944 ph25

p6*4.111 ph23

lo to 4 times l1

;end DIPSI2

p16:gp3

d16

(p21:sp21 ph1):f1

DELTA

(center (p22:sp22 ph1):f1 (p22:sp23 ph6):f2)

DELTA1 pl1:f1

(ralign (p1 ph1):f1 (p21:sp24 ph9):f2)

DELTA2

d0

(p22:sp22 ph5):f1

d0

DELTA2 pl1:f1

(lalign (p1 ph7):f1 (p25:sp25 ph10):f2)

DELTA3

(p22:sp22 ph1):f1

3u

p16:gp1*EA1

d16

(p21:sp26 ph10):f2

d4

(center (p22:sp22 ph1):f1 (p22:sp23 ph6):f2)

3u

p16:gp1*EA2

d16 pl12:f2

DELTA4 BLKGRAD

go=2 ph31 cpd2:f2

d1 do:f2 mc #0 to 2

F1EA(igrad EA1 & igrad EA2, id0 & ip9*2 & ip6*2 & ip31*2)

exit

ph1=0

ph2=1

ph5=0 0 2 2

ph6=0

ph7=3

ph8=3

ph9=2 0

ph10=0 0 0 0 2 2 2 2

ph23=3

ph25=1

ph31=0 2 0 2 2 0 2 0

;pl0 : 0W

;pl1 : f1 channel - power level for pulse (default)

;pl2 : f2 channel - power level for pulse (default)

;pl10 : f1 channel - power level for mixing sequence (DIPSI-2)

;pl12: f2 channel - power level for CPD/BB decoupling

;p1 : f1 channel - 90 degree high power pulse

;p2 : f1 channel - 180 degree high power pulse

;p6 : f1 channel - 90 degree pulse for mixing sequence (DIPSI-2)

;p21 : 550u excitation pulses

;p22 : 600u UR and Inversion pulses

;p25 : 1100u 13C refocusing pulse

;p16: homospoil/gradient pulse

;spnam21 : jc01_BEBOP_zx_550u_BW10_RF20_pm20_Hc0.99997119.pul

;spnam22 : jc02_BURBOP_x_600u_BW10_RF20_pm20_matched.pul

;spnam23 : jc03_BIBOP_600u_BW37.5_RF10_pm5_matched.pul

;spnam24 : jc05_BEBOP_zy_550u_BW37.5_RF10_pm5_matched.pul

;spnam25 : jc07_BURBOP_y_1100u_BW37.5_RF10_pm5_Hc0.999876221.pul

;spnam26 : jc09_BEBOP_-yz_550u_BW37.5_RF10_pm5_matched.pul

;sp21 : 20 kHz Rf Amplitude

;sp22 : 20 kHz Rf Amplitude

;sp23 : 10 kHz Rf Amplitude

;sp24 : 10 kHz Rf Amplitude

;sp25 : 10 kHz Rf Amplitude

;sp26 : 10 kHz Rf Amplitude

;d0 : incremented delay (2D) [3 usec]

;d1 : relaxation delay; 1-5 * T1

;d4 : 1/(4J)XH

;d3 : reduced d4

;d9 : TOCSY mixing time

;d11: delay for disk I/O [30 msec]

;d16: delay for homospoil/gradient recovery [200ms]

;d30: delay before experiment starts (use if running extended multizgs)

;cnst4: = J(XH)

;cnst3: reduced cnst4 (’ernst angle’)

;inf1: 1/SW(X) = 2 * DW(X)

;in0: 1/(2 * SW(X)) = DW(X)

;nd0: 2

;NS: 1 * n

;DS: >= 16

;td1: number of experiments

;FnMODE: echo-antiecho

;cpd2: decoupling according to sequence defined by cpdprg2

;pcpd2: f2 channel - 90 degree pulse for decoupling sequence

;use gradient ratio: gp 1 : gp 2 : gp 3

; 80 : 33 : 43

;for z-only gradients:

;gpz1: 80%

;gpz2: 33% ; gradient before mixing

;gpz3: 43% ; gradient after mixing

;use gradient files:

;gpnam1: SMSQ10.100

;gpnam2: SMSQ10.100

;gpnam3: SMSQ10.100

##  S9 EPSI

All spectra in this section were acquired on a $600\text{\,}\mathrm{MHz}$
Bruker AVIII spectrometer, equipped with a Prodigy N2 broadband cryoprobe with
a nominal maximum $z$-gradient strength of
$66\text{\,}\mathrm{G}\text{\,}{\mathrm{cm}}^{-1}$. The sample used was
$45\text{\,}\mathrm{m\textsc{m}}$ ferulic acid in DMSO-$d_{6}$.

###  S9.1 Cost function

The cost function used for EPSI optimisations (called epsi_gradient_drift)
comes pre-installed with POISE. It is more involved than usual, as the
entirety of the EPSI processing must be carried out within the cost function
(POISE itself only providing a function to get the FID). Briefly, this
processing consists of the following steps:

  * •

Removal of the “group delay” (by circularly shifting the FID).

  * •

Reshaping of the 1D FID into a 2D $(k,t_{2})$ matrix.

  * •

Removing parts of the FID acquired under negative gradients, as well as during
the delay between gradients.

  * •

Absolute-value processing.

  * •

Apodisation.

  * •

Removing any rows (i.e. any values of $t_{2}$) for which the maximum intensity
is below a certain threshold. This ensures that only parts of the $(k,t_{2})$
matrix with strong signal are used for the calculation of the slope.

The slope of $k$ versus $t_{2}$ can then be calculated by finding, for each
row, the value of $k$ where the intensity is maximum, and performing a least-
squares fitting. Technically $k/k_{\text{max}}$, not $k$, is used instead of
$k$, but $k_{\text{max}}$ is merely a constant factor which has no influence
on the results.

For more details, please refer to the source code which is attached below:

⬇

”””

costfunctions.py

—————-

Contains default cost functions. Custom cost functions should be added in

costfunctions_user.py (using the same format).

You can add custom cost functions in here (they will work), but they will be

shadowed by any similarly named cost functions in costfunctions_user.py. Also,

this file will be overwritten if POISE is reinstalled.

SPDX-License-Identifier: GPL-3.0-or-later

”””

import numpy as np

from .cfhelpers import *

from .shared import _g

def noe_1d():

”””

Measures the intensity of peaks in the spectrum, *except* for anything

within 25 Hz of the selectively excited peak.

”””

# SPOFFS2 plus O1 is the frequency of the selective pulse, in Hz.

# We can get the values of both using getpar(), which automatically

# converts numeric parameters to floats.

f = getpar(”SPOFFS2”) + getpar(”O1”)

# Define the (rough) bandwidth of the selective pulse. We will ignore the

# region of the spectrum that is bw Hz wide and centred on the selective

# pulse.

bw = 50.0

# To get the spectrum *without* the ignored region, we split up the

# spectrum into two ”portions”. We need to construct the ’bounds’ parameter

# for both portions.

sfo1 = getpar(”SFO1”)

lowerhalf = f”{(f + bw/2)/sfo1:.3f}..” # frequencies (f + bw)/2 and above

upperhalf = f”..{(f - bw/2)/sfo1:.3f}” # frequencies up to (f - bw)/2

# We use get1d_real() to get both portions of the spectrum, passing the

# bounds as appropriate. The intensities of the spectrum are obtained by

# summing all points (np.sum).

# Note that manually passing the bounds parameter will override the F1P/F2P

# parameters if set by the user (this is probably a good thing).

upper_integral = np.sum(get1d_real(bounds=upperhalf))

lower_integral = np.sum(get1d_real(bounds=lowerhalf))

# Lastly, we take the negative absolute value of the resulting sum, in

# order to make sure that the resulting cost function is negative

# regardless of how the NOE peaks are phased (positive or negative).

return -abs(upper_integral + lower_integral)

def minabsint():

”””

Cost function which minimises the absolute (magnitude-mode) intensity of

the spectrum. This is probably the easiest cost function. :-)

”””

# Get the real part of the spectrum as a numpy ndarray. Because no bounds

# are explicitly passed, this defaults to the whole spectrum if F1P/F2P are

# not specified. If F1P/F2P *are* specified (either via the ‘dpl‘ TopSpin

# command, or manually), this will only return the part of the spectrum

# between those two values.

r = get1d_real()

# …and the imaginary part

i = get1d_imag()

# The magnitude-mode spectrum is obtained by combining these.

mag = np.abs(r + 1j * i)

# The intensity of the magnitude-mode spectrum is just the sum of all

# points.

return np.sum(mag)

def maxabsint():

”””

Cost function which maximises the absolute (magnitude-mode) intensity of

the spectrum.

”””

r = get1d_real()

i = get1d_imag()

# This is the same as minabsint except that we have a negative sign.

# Because the optimisation always seeks to *minimise* the cost function,

# this essentially tries to *maximise* np.sum(…), i.e. maximise the

# spectral intensity.

return -np.sum(np.abs(r + 1j*i))

def minrealint():

”””

Minimises the intensity of the real part of the spectrum. If the spectrum

has negative peaks this cost function will try to maximise those.

”””

return np.sum(get1d_real())

def maxrealint():

”””

Maximises the intensity of the real part of the spectrum.

”””

return -np.sum(get1d_real())

def zeronetrealint():

”””

Tries to get the *net* intensity of the real spectrum to be as close to

zero as possible. Note that dispersion-mode peaks will contribute to this

cost function (as they add to zero).

”””

return np.abs(np.sum(get1d_real()))

def zerorealint():

”””

Tries to get the real spectrum to be as close to zero as possible. This

takes the absolute value of the real spectrum before summing, so

dispersion-mode peaks do not add to zero.

”””

return np.sum(np.abs(get1d_real()))

def zerorealint_squared():

”””

Tries to get the real spectrum to be as close to zero as possible. This

uses a sum-of-squares metric, which (compared to “zerorealint“) leads to

two main differences. Firstly, dispersion-mode peaks do not sum to zero

because sign information is lost upon squaring. Secondly, it more strongly

penalises large deviations away from zero.

”””

s = get1d_real()

return np.sum(s * s)

def epsi_gradient_drift():

”””

Calculates the amount of ’gradient drift’ seen in a 1D EPSI acquisition, as

reflected by the position of the ’echo’ moving over time. This can be

caused by imbalanced positive and negative gradients.

”””

from numpy.polynomial.polynomial import Polynomial

# — Calculate key parameters ——————————-

fid = get1d_fid(remove_grpdly=True)

# Number of complex points in k-space, i.e. number of points per EPSI

# gradient. TD2 is total number of real & imag points in FID. L3 is number

# of EPSI loops (one loop includes both pos + neg gradient). So TD2 / (2 *

# L3) is the number of points in one EPSI gradient (i.e. only positive

# gradient). The extra factor of 2 is because we’re interested only in

# complex points.

td_k = int(getpar(”TD”) / (2 * 2 * getpar(”L3”)))

# Time between consecutive EPSI positive gradients, in seconds

td_t2 = int(getpar(”L3”))

dw_eff = getpar(”AQ”) / td_t2

# — Perform EPSI processing ——————————–

# Reshape into 2D matrix

ser = fid.reshape((-1, td_k))

# Discard the part acquired with negative gradients

ser = ser[0::2, :]

# Discard any part of the spectrum that was not acquired during an EPSI

# gradient, i.e. if the delay D6 was nonzero.

if getpar(”D6”) > 0:

ineligible_epsi_points = int(1e6 * getpar(”D6”) / (getpar(”DW”) * 2))

td_k = td_k - ineligible_epsi_points

ser = ser[:, :td_k]

# — Apodisation ——————————————–

# Along k-dimension (Hamming window)

alpha_0 = 0.54

ks = np.linspace(0, 1, td_k)

k_winfunc = (alpha_0

- (1 - alpha_0) * np.cos(2 * np.pi * np.linspace(0, 1, td_k)))

ser = ser * k_winfunc[np.newaxis, :]

# Along direct dimension

t2_winfunc = np.sin(np.pi * np.linspace(0, 1, td_t2))

ser = ser * t2_winfunc[:, np.newaxis]

abs_ser = np.abs(ser)

# Calculate k- and t2-axes

t2_values = np.arange(td_t2) * dw_eff

k_values = np.linspace(-0.5, 0.5, td_k)

# Drop all rows (i.e. all values of t2) for which the maximum is less than

# 20% of the overall maximum.

threshold_amp = 0.2 * np.max(abs_ser)

maxima_along_rows = np.max(abs_ser, axis=1)

indices_to_use = np.nonzero(maxima_along_rows > threshold_amp)

# Locate the maxima

maximal_indices_along_k = np.argmax(abs_ser, axis=1)

# Get the value of k at each maximum

k_max_values = k_values[maximal_indices_along_k]

poly = Polynomial.fit(x=t2_values[indices_to_use],

y=k_max_values[indices_to_use], deg=1)

return np.abs(poly.coef[1])

# The cost functions below were used in the POISE paper, but are less likely
to

# be generally useful. If you want to use them, simply uncomment them (remove

# the three single quotes around them).

”’

def asaphsqc():

”””

Maximises the intensity of the projection of a 2D spectrum onto the f2

axis.

”””

# First we get the ’rr’ part of the spectrum (the displayed part). This

# assumes that it has been phased already (by the AU program).

# This is a 2D numpy array:

# [ [1 , 2, 3, …, si2],

# [. , ., ., …, …],

# [. , ., ., …, …],

# [si1, ., ., …, …]]

# where si1 and si2 are the values of SI in the indirect and direct

# dimension respectively.

spec = get2d_rr()

# We take the ’skyline’ projection, i.e. the highest point along each

# column. This is basically equivalent to TopSpin’s ‘f2projp‘ command.

# If you want a projection onto f1 instead of f2, you can use axis=1.

proj = np.amax(spec, axis=0)

return -np.sum(proj)

”’

”’

def specdiff():

”””

Returns the ”spectral difference” between the current spectrum and a

reference spectrum with EXPNO 1 and PROCNO 1\.

”””

# We can construct the path to the reference spectrum using make_p_spec().

reference_path = make_p_spec(expno=1, procno=1)

# The way this is currently coded, the reference spectrum has to be placed

# in EXPNO 1 and PROCNO 1. In theory, this can be specified using a TopSpin

# parameter in the optimisation dataset, for example ”cnst30” (or any other

# random cnst). In order to use this, just add:

# expno = # int(getpar(”cnst30”))

# and then after that you can use

# make_p_spec(expno=expno, procno=1)

# Then we can get the reference, or ’target’, spectrum as a numpy array.

# Note that this will respect the F1P/F2P parameters, *not* of the

# reference spectrum, but of the spectrum being optimised. This is because

# get1d_real() does not actually read F1P and F2P from the spectrum. It

# just uses the global values stored in _g.spec_f1p and _g.spec_f2p, which

# refer to the F1P/F2P values of the spectrum being optimised.

target = get1d_real(p_spec=reference_path)

# … and the current spectrum that is being optimised.

spec = get1d_real()

# Finally, return the 2-norm of the difference between the two normalised

# spectra.

return np.linalg.norm(target/np.linalg.norm(target) \-

spec/np.linalg.norm(spec))

”’

”’

def dosy():

”””

Tries to make a spectrum (1D diffusion experiment with variable gradient

strength) have 25% intensity of a target spectrum (1D diffusion experiment

with minimum gradient strength). This should be used in conjunction with

the dosy_opt Python programme (after POISE is installed via ‘pip‘, this can

be installed from the TopSpin command-line via ‘poise –install dosy‘.)

”””

# The target spectrum is hardcoded to have EXPNO 99998. See psyche() for

# more discussion about this.

# As in psyche(), the F1P/F2P parameters are not taken from the reference

# spectrum but from the spectrum being optimised.

reference_path = make_p_spec(expno=99998, procno=1)

target = get1d_real(p_spec=reference_path)

# This is the spectrum being optimised.

spec = get1d_real()

# The intensity of the optimised spectrum is np.sum(spec), and likewise for

# the target spectrum. We take the ratio of the two and measure how ”far

# away” this is from the ideal value of 0.25 (which corresponds to 75%

# attenuation).

return np.abs(np.sum(spec)/np.sum(target) - 0.25)

”’

”’

def dosy_aux():

”””

Non-absolute value of dosy(). To be used in the first stage of dosy_opt.

”””

reference_path = make_p_spec(expno=99998, procno=1)

target = get1d_real(p_spec=reference_path)

spec = get1d_real()

return np.sum(spec)/np.sum(target) - 0.25

”’

”’

def dosy_2p():

”””

Two-parameter DOSY optimisation.

”””

reference_path = make_p_spec(expno=99998, procno=1)

target = get1d_real(p_spec=reference_path)

spec = get1d_real()

return np.abs(np.sum(spec)/np.sum(target) - 0.25) + getpar(”D20”)

”’

###  S9.2 Reference grid search

The parameter optimised here is cnst16, which represents the ratio of nominal
negative gradient amplitude to nominal positive gradient amplitude during the
EPSI acquisition. In the main text, this is referred to as $\alpha$. On an
ideal instrument, the optimised value of cnst16 would simply be 1.

Figure S8:  Dependence of the cost function on the value of cnst16. The
minimum occurs at around cnst16 = 1.0004.

###  S9.3 Optimisation

Entry | Algorithm | Optimum found | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---
1 | Nelder–Mead | 1.000375 | 10 | 49–52
2 | MDS | 1.000375 | 10 | 49–50
3 | BOBYQA | 1.000383–1.000388 | 7 | 34–35
Table S13:  Results of EPSI gradient imbalance optimisations. The POISE
routine used here is: . The pulse programme used was zg_epsi.2_jy (Section
S9.4). Key acquisition parameters: NS=1, DS=0, L4=0, D1=0.5, DW=0.7 (note
units of DW are $\mathrm{\SIUnitSymbolMicro s}$), AQ=0.092, GPZ15=35.644 (the
maximum gradient amplitude on this machine is
$67.5\text{\,}\mathrm{G}\,\mathrm{c}\mathrm{m}^{-1}$). To avoid hardware
damage AQ should always be kept short, ideally under
$100\text{\,}\mathrm{ms}$. Note that such a short value of D1 is only possible
because there is always an extra delay between successive experiments required
for pulse programme compilation, etc. A short recovery delay should not be
used when DS > 0 or NS > 1. $n_{\text{fev}}$ refers to the number of function
evaluations, i.e. the number of spectra acquired in the course of the
optimisation. The results are aggregated from five separate runs for each
optimisation algorithm.

###  S9.4 Excitation–EPSI pulse programme

⬇

;zg_epsi.2_jy

;

;download from: https://git.io/JYGEa

;

;90 pulse followed by EPSI acquisition, with prephasing gradient

;modified from uftocsy, Giraudeau et al.

;https://madoc.univ-nantes.fr/course/view.php?id=24710&section=1

;IMPORTANT: set p15 + d6 = DW x TD(F2)/(2xL3)

;$CLASS=HighRes

;$DIM=1D

;$TYPE=

;$SUBTYPE=

;$COMMENT=

#include <Avance.incl>

#include <Grad.incl>

#include <De.incl>

#include <Delay.incl>

”d0=3u”

”d11=30m”

”p15=(td*dw/(2*l3))-d6”

”cnst17=cnst16” ; so that cnst16 shows up in ased

1 ze

”cnst17=cnst17” ; so that cnst16 shows up in ased

d11

30m zd

; dummy scans

4 30m

10u

4u

d1

4u

100u UNBLKGRAD

20u pl1:f1

(p1 ph1):f1

10u

p15:gp25

10u

5 p15:gp15

d6

p15:gp16

d6

lo to 5 times l3

p27:gp27

100u BLKGRAD

lo to 4 times l4

; end of dummy scans

2 30m

10u

4u

d1

4u

100u UNBLKGRAD

20u pl1:f1

; excitation

(p1 ph1):f1

; prephasing gradient

10u

p15:gp25

10u

; EPSI acquisition

ACQ_START(ph30,ph31)

1u DWELL_GEN:f1

3 p15:gp15

d6

p15:gp16*cnst16

d6

lo to 3 times l3

100u BLKGRAD

rcyc=2

100u UNBLKGRAD

p27:gp27

100u BLKGRAD

30m mc #0 to 2 F0(zd)

d17

exit

ph1=0

ph2=0 1

ph3=0

ph29=0

ph30=0

ph31=0 2

;pl1: f1 channel - power level for pulse (default)

;p1: f1 channel - 90 degree high power pulse

;d1: relaxation delay; 1-5 * T1

;d6: delay between EPSI gradients

;d17: delay after experiment (not necessary)

;p15: EPSI gradient duration

;gpz15: EPSI positive gradient

;gpz16: EPSI negative gradient

;gpz25: prephasing gradient

;gpz27: purge gradient

;ns: 1

;l3: number of loops for EPSI acquisition

##  S10 PSYCHE

All spectra in this section were acquired on a $600\text{\,}\mathrm{MHz}$
Bruker AVIII spectrometer, equipped with a Prodigy N2 broadband cryoprobe with
a nominal maximum $z$-gradient strength of
$66\text{\,}\mathrm{G}\text{\,}{\mathrm{cm}}^{-1}$. The sample used was
$45\text{\,}\mathrm{m\textsc{m}}$ andrographolide in DMSO-$d_{6}$.

###  S10.1 Cost function

The cost function, named specdiff, is given by:

$f=\left|\left|\frac{\mathbf{S}}{||\mathbf{S}||}-\frac{\mathbf{T}}{||\mathbf{T}||}\right|\right|$
(5)

where $\mathbf{S}$ and $\mathbf{T}$ are vectors corresponding to the spectrum
under optimisation and a target spectrum respectively. The target spectrum was
acquired using a proton pulse-acquire sequence (zg). The optimisation spectrum
was acquired using a $None$–delay–refocusing element–delay–acquire sequence,
as described in the main text. (This pulse programme is attached further
below.) The spectral region used in the optimisation was restricted to
$1.15\text{\,}\mathrm{ppm}$ and above, in order to exclude methyl singlets.
This is not strictly necessary (roughly the same results will be obtained if
the entire spectral window is used), but helps to reduce noise in the cost
function.

specdiff is not pre-installed with POISE: to enable it, uncomment the
specdiff() function in poise_backend/costfunctions.py.

###  S10.2 Flip angle reference grid search

Due to the large amount of time required for multiple-parameter grid searches,
we chose to only run this over a single parameter, namely the flip angle
$\beta$.

Figure S9: Dependence of cost function on PSYCHE flip angle (cnst20).

The cost function is minimised at around 15–$None$. We can use this
opportunity to have a closer look at exactly how this cost function works. It
penalises the low signal-to-noise seen at smaller flip angles: in this case
$||\mathbf{S}||$ is small and consequently the noise in
$\mathbf{S}/||\mathbf{S}||$ has a larger intensity. This contributes to a
larger value of $\mathbf{S}/||\mathbf{S}||-\mathbf{T}/||\mathbf{T}||$ and
hence a larger value of the cost function. On the other hand, at larger flip
angles this is not a problem: instead one starts to observe phase distortions
in the peaks. This causes
$\mathbf{S}/||\mathbf{S}||-\mathbf{T}/||\mathbf{T}||$ to be larger in the
vicinity of the peaks, which also pushes $f$ up. The cost function is
therefore minimised in a “compromise” region which has acceptable signal-to-
noise but also not too much phase distortion. This is depicted in Figure S10.

Figure S10: A closer look at the cost function at different flip angles. Left
column: The target spectrum $\mathbf{T}$ (this is the same for all three
rows). Middle column: The 1D test spectrum $\mathbf{S}$, acquired using
different flip angles (the optimum is around $None$). Right column: The vector
$\mathbf{S}/||\mathbf{S}||-\mathbf{T}/||\mathbf{T}||$. One can see how at low
flip angles, the main contribution to the cost function is noise. At higher
flip angles, the noise contributes less and the increases in cost function are
due to phase distortion.

###  S10.3 Optimisations

In this section we present the overall results of all the PSYCHE optimisations
ran (with 1, 2, 3, and 4 parameters being optimised at once). From the
aggregated results (Table S14), it is clear that increasing the number of
parameters leads to larger decreases in the cost function, i.e. better
performance as a refocusing element, although there is a drawback in that more
time is required for convergence. We note here that an increase in the number
of parameters leads to greater complexity in the cost function. This increases
the likelihood of converging to a local minimum which may be suboptimal. The
tradeoff between various factors such as potentially improved performance, the
time required for optimisation, and the chances of being caught in a shallow
local minimum must be carefully weighed when deciding how complicated an
optimisation to run. In general, our opinion is that optimising two parameters
(the flip angle and gradient amplitude) provides a good balance between all of
these competing factors, and is most suitable for routine use.

First, we present an overview of all the optimisations run.

| Aggregated results from all runs | Optimised parameters from best run
---|---|---
Description | $f_{\mathrm{best}}$ | $n_{\mathrm{fev}}$ | Time (s) |  $\beta$ (∘) |  $g$ (%) |  $\Delta F$ (kHz) |  $\tau_{\mathrm{p}}$ (ms)
Initial point | 0.353 | – | – | (25.0) | (2.00) | (10.0) | (30.0)
1 parameter | 0.340–0.343 | 5–10 | 84–168 | 17.5 | (2.00) | (10.0) | (30.0)
2 parameters | 0.325–0.338 | 12–33 | 205–565 | 16.5 | 1.00 | (10.0) | (30.0)
3 parameters | 0.320–0.328 | 18–77 | 315–1344 | 17.4 | 1.73 | 15.9 | (30.0)
4 parameters | 0.316–0.333 | 39–85 | 705–1504 | 13.9 | 1.45 | 15.9 | 36.0
Table S14:  Summary of optimisations on PSYCHE refocusing element. Details of
the routines used are described in Tables S15, S16, S17 and S18. Parameters in
parentheses were not optimised (they are inherited from the initial point).
$n_{\text{fev}}$ refers to the number of function evaluations, i.e. the number
of spectra acquired in the course of the optimisation. The results are
aggregated from five separate runs for each optimisation algorithm. The
optimised parameters in each row are taken from the best run, as judged by the
lowest value of the cost function found.

We follow this with details of the individual classes of optimisation. For the
1-parameter optimisation, the results are summarised in a similar format to
previous tables. However, for the remainder of the optimisations, we present
only summaries of $n_{\text{fev}}$ and the time taken. For details of the
optima found, we point the interested reader in the direction of the poise.log
files, as well as the parse_log() function to help with parsing these (please
refer to the documentation for more details).

Entry | Algorithm | Optimum found ($\mathrm{\SIUnitSymbolDegree}$) | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---
1 | Nelder–Mead | 15.0–20.0 | 8–10 | 135–168
2 | MDS | 17.5–20.0 | 8 | 133–136
3 | BOBYQA | 18.4–19.9 | 5 | 84–85
Table S15:  PSYCHE 1-parameter (flip angle) optimisations. The POISE routine used was: . Key acquisition parameters: NS=2, DS=1, D1=1.5.  Entry | Algorithm | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---
1 | Nelder–Mead | 18–24 | 307–410
2 | MDS | 25–33 | 424–565
3 | BOBYQA | 12–16 | 205–271
Table S16:  PSYCHE 2-parameter (flip angle and gradient amplitude) optimisations. The POISE routine used was: . Key acquisition parameters: NS=2, DS=1, D1=1.5.  Entry | Algorithm | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---
1 | Nelder–Mead | 33–43 | 576–770
2 | MDS | 46–77 | 804–1344
3 | BOBYQA | 18–31 | 315–540
Table S17:  PSYCHE 3-parameter (flip angle, gradient amplitude, and bandwidth) optimisations. The POISE routine used was: . Key acquisition parameters: NS=2, DS=1, D1=1.5.  Entry | Algorithm | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---
1 | Nelder–Mead | 40–47 | 733–845
2 | MDS | 57–85 | 1006–1504
3 | BOBYQA | 39–62 | 705–1130
Table S18:  PSYCHE 4-parameter (flip angle, gradient amplitude, bandwidth, and
pulse duration) optimisations. The POISE routine used was: . Key acquisition
parameters: NS=2, DS=1, D1=1.5.

The best of each of these optimisations are shown in the next figure.

Figure S11:  Insets of PSYCHE spectra of a sample of
$45\text{\,}\mathrm{m\textsc{m}}$ andrographolide in DMSO-$d_{6}$, obtained
before and after optimisation of up to 4 parameters. The original coupled ^1H
spectrum is shown as a reference. The time required for each optimisation is
indicated on each spectrum.

###  S10.4 AU programme used for optimisations

The poise_psyche AU programme was used for all the optimisations above. It is
the same as poise_1d, except that it calls on a separate Python script to
create a new double saltire pulse with the specified parameters prior to
acquisition (similar to TopSpin’s WaveMaker software):

⬇

XCMD(”sendgui xpy make_double_saltire”)

ZG

EFP

APBK

QUIT

The Python script make_double_saltire is described next.

###  S10.5 Python script for generating double saltire

Please note that this is a Python script to be run inside TopSpin. It can be
downloaded from the link in the docstring.

⬇

”””

make_double_saltire.py

———————-

Download from: https://git.io/JUHOZ

This script reads in a pulse duration (stored in us as P40) and bandwidth

(stored in Hz as CNST21), then creates a double saltire shaped file called

”auto_double_saltire” which can be directly used as the shaped file in

PSYCHE experiments.

This can be included in an AU programme using the line:

XCMD(”sendgui xpy make_double_saltire”)

SPDX-License-Identifier: GPL-3.0-or-later

”””

from __future__ import division

from math import sin, cos, pi, sqrt, degrees, atan2

def make_double_saltire(duration, bandwidth):

”””

Parameters

———-

duration : float

Time in units of seconds.

bandwidth : float

Bandwidth in units of Hz.

Returns

——-

A : list of float

Amplitudes, scaled between 0 and 100\.

phi : list of float

Phases, in degrees between -180 and 180\.

These lists can be directly passed to the Bruker SAVE_SHAPE routine, as

demonstrated in this script. The amplitude of the pulse should be

automatically calculated in the pulse programme.

”””

# Our strategy is to generate the first saltire first.

# Number of points in half of the pulse. We use 1 point per us.

npoints = int(duration * 1000000 / 2)

# 20% smoothing.

s = 0.2

# This is equivalent to ls = np.linspace(0, 1, npoints), but you can’t use

# numpy in TopSpin…

ts = [i/(npoints - 1) for i in range(0, npoints)]

# Construct the amplitude piecewise function for a single chirp.

A_chirp = [sin((pi * t) / (2 * s)) if t < s

else 1 if s <= t and t <= 1 - s

else sin((pi * (1 - t)) / (2 * s))

for t in ts]

# These are phases for a single chirp, not a saltire.

phi_chirp = [pi * (duration/2) * bandwidth * ((t - 0.5) ** 2) for t in ts]

# Calculate x-coefficients (Cx) for the saltire.

Cx = [100 * a * cos(phi) for a, phi in zip(A_chirp, phi_chirp)]

# The Cy’s are 100 * a * sin(phi) …, but that’s only for one chirp. For a

# saltire, the two component chirps have equal and opposite Cy’s (because

# sin(phi) = -sin(-phi)), so the total Cy is simply zero for every point.

# In principle we could make a list full of zeroes, like Cy = [0 for _ in

# range(npoints)], but that’s a waste of memory and time.

# Incidentally, the Cx’s for both chirps add up since cos(phi) = cos(-phi),

# but because the amplitudes are scaled to [0, 100], there’s no point

# actually adding them up.

# Convert back to polar coordinates. Ordinarily A = sqrt(xi * xi + yi *

# yi), but since all the yi’s are zero, we just leave it out of the

# equation.

A = [sqrt(xi * xi) for xi in Cx]

# Likewise here we use atan2(0, xi) instead of atan2(yi, xi).

phi = [degrees(atan2(0, xi)) for xi in Cx]

# Finally, we double the lists (because it’s a double saltire).

return A * 2, phi * 2

if __name__ == ”__main__”:

# Let the user know what’s going on

SHOW_STATUS(”make_double_saltire: generating saltire…”)

# Pulses are by default in us, so must convert first

duration = float(GETPAR(”P 40”))/1000000

# Bandwidth is specified in Hz.

bandwidth = float(GETPAR(”CNST 21”))

# Generate the amplitudes and phases.

A, phi = make_double_saltire(duration, bandwidth)

# Save it to disk.

shape_name = ”auto_double_saltire”

SAVE_SHAPE(shape_name, ”Excitation”, A, phi)

# Use it as SPNAM40. It’s too easy to forget this!

PUTPAR(”SPNAM 40”, shape_name)

###  S10.6 1D spin echo pulse programme

⬇

;psyche_1dopt_jy

;

;download from: https://git.io/JUHOS

;

;90 - 15m delay - double saltire - 15m delay - acquire

;Jonathan Yong - University of Oxford - 26 Aug 2020

;Optimisation parameters:

;cnst20 : flip angle (degrees)

;gpz10 : gradient amplitude (%)

;cnst21 : bandwidth (Hz)

;p40 : saltire + gradient duration (us)

#include <Avance.incl>

#include <Grad.incl>

”d2=15m”

”cnst10=(cnst20/360)*sqrt((2*cnst21)/(p40/2000000))” ; rf amplitude of saltire

”cnst11=1000000/(p1*4)” ; rf amplitude hard pulse

”spw40=plw1*(cnst10/cnst11)*(cnst10/cnst11)” ; power level of saltire

”acqt0=0”

baseopt_echo

1 ze

2 30m

50u BLKGRAD

d1 pl1:f1

50u UNBLKGRAD

p1 ph1

d2

d16

p16:gp1

d16

(center (p40:sp40 ph3):f1 (p40:gp10) )

d16

p16:gp1

d16

d2

go=2 ph31

30m mc #0 to 2 F0(zd)

50u BLKGRAD

exit

ph1 = 0

ph3 = 0 1 2 3

ph31= 0 2 0 2

;pl1 : f1 channel - power level for pulse (default)

;p1 : f1 channel - high power pulse

;p40: shaped pulse in PS element

;p21: gradient during PS element

;d1 : relaxation delay; 1-5 * T1

;d2 : delay on either side of double saltire

;d16: gradient recovery delay

;p16: CTP gradient duration

;ns: 1 * n, total number of scans: NS * TD0

;cnst20: flip angle (degrees)

;cnst21: bandwidth of PSYCHE element (Hz) (ca. 10000)

;gpnam10: RECT.1

;gpz10 : 2%

;gpnam1 : SINE.100

;gpz1 : 50%

##  S11 1D NOESY-based solvent suppression

All spectra in this section were acquired on a $400\text{\,}\mathrm{MHz}$
Bruker Avance NEO spectrometer, equipped with a broadband Smart probe with a
nominal maximum $z$-gradient strength of
$50\text{\,}\mathrm{G}\text{\,}{\mathrm{cm}}^{-1}$. A sample of rodent urine
in D2O was used, kindly provided by Fay Probert and Abi Yates (Department of
Chemistry, University of Oxford).

###  S11.1 Cost function

The cost function used for all optimisations here was zerorealint_squared,
which attempts to minimise the squared intensity of the real spectrum:

$f=\sum_{i}\mathbf{S}_{\mathrm{real},i}^{2}$ (6)

The aim of this is to strongly penalise large deviations from the baseline,
which are more likely to potentially obscure nearby peaks. This cost function
comes pre-installed with POISE. The region of the spectrum under consideration
was restricted to between 4.65 and $4.75\text{\,}\mathrm{ppm}$ using the dpl
function in TopSpin.

###  S11.2 Transmitter offset reference grid search

It is not feasible to run a grid search over four parameters at the same time.
As with the PSYCHE section, we therefore only ran one grid search for the
first parameter, i.e. the transmitter offset O1.

Figure S12:  Insets of the residual water peak, and the value of the cost
function thus calculated, as a function of transmitter offset (O1).

It is evident here that there are multiple possible minima which POISE could
optimise to. POISE is not a global minimiser; it will only converge to a local
minimum. In fact, a similar situation was observed when we ran grid searches
over the other three parameters: multiple potential minima were observed, due
to an oscillatory dependence on the parameters under optimisation.

###  S11.3 Optimisations

It is likely that there are multiple different sets of optimised parameters
which give similar minima in the cost function. We assume that the user does
not mind exactly which minimum is obtained, as long as there is a sufficient
decrease in the cost function. In other words, we assume that the user is not
concerned with exactly which combination of parameters provides the best
solvent suppression (as long as they are within reasonable bounds, e.g.
presaturation power is not too high; these should be set appropriately in the
POISE routine).

As with the PSYCHE optimisations (Section S10), we note that optimising more
parameters simultaneously takes a longer time and risks being caught in local
minima. However, it is also clear that this has the potential of leading to a
greater decrease in the cost function. In our experience, optimising two
parameters (O1 and CNST20) is usually sufficient to give a good decrease in
the cost function (i.e. good water suppression). For routine usage this two-
parameter optimisation should be considered.

The detailed optimisation results are presented below, in a similar format to
the PSYCHE section (Section S10.3). On the $400\text{\,}\mathrm{MHz}$
spectrometer we used, the default O1 value of $1880.61\text{\,}\mathrm{Hz}$
corresponds to exactly $4.70\text{\,}\mathrm{ppm}$.

| Aggregated results from all runs | Optimised parameters from best run
---|---|---
Description | $f_{\mathrm{best}}/10^{18}$ | $n_{\mathrm{fev}}$ | Time (s) |  O1 ($\mathrm{Hz}$) |  CNST20 ($\mathrm{Hz}$) |  D8 ($\mathrm{s}$) |  D1 ($\mathrm{s}$)
Initial point | 14.7 | – | – | (1880.61) | (50.0) | (0.100) | (2.00)
1 parameter | 1.85–2.49 | 6–12 | 259–520 | 1880.41 | (50.0) | (0.100) | (2.00)
2 parameters | 1.21–9.68 | 9–21 | 390–911 | 1880.20 | 51.94 | (0.100) | (2.00)
3 parameters | 1.24–4.20 | 19–26 | 825–1128 | 1879.90 | 47.79 | 0.118 | (2.00)
4 parameters | 0.165–1.65 | 25–53 | 1143–2314 | 1881.10 | 53.28 | 0.150 | (3.00)
Table S19:  Summary of optimisations on 1D NOESY / presaturation pulse sequence for solvent suppression.. Details of the routines used are described in Tables S20, S21, S22 and S23. Parameters in parentheses were not optimised (they are inherited from the initial point). $n_{\text{fev}}$ refers to the number of function evaluations, i.e. the number of spectra acquired in the course of the optimisation. The results are aggregated from five separate runs for each optimisation algorithm. The optimised parameters in each row are taken from the best run, as judged by the lowest value of the cost function found.  Entry | Algorithm | Optimum found ($\mathrm{Hz}$) | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---|---
1 | Nelder–Mead | 1880.24–1880.49 | 10–12 | 436–519
2 | MDS | 1880.24–1880.36 | 12 | 518–520
3 | BOBYQA | 1880.34–1880.47 | 6–7 | 259–303
Table S20:  Solvent suppression 1-parameter (transmitter offset) optimisations. The POISE routine used was: . Key acquisition parameters: NS=2, DS=4, AQ=4.19.  Entry | Algorithm | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---
1 | Nelder–Mead | 17–21 | 737–911
2 | MDS | 15 | 649–652
3 | BOBYQA | 9–14 | 390–608
Table S21:  Solvent suppression 2-parameter (transmitter offset and presaturation power) optimisations. The POISE routine used was: . Key acquisition parameters: NS=2, DS=4, AQ=4.19.  Entry | Algorithm | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---
1 | Nelder–Mead | 23–26 | 1002–1128
2 | MDS | 24–25 | 1034–1080
3 | BOBYQA | 19–26 | 825–1133
Table S22:  Solvent suppression 3-parameter (transmitter offset, presaturation power, and mixing time) optimisations. The POISE routine used was: . Key acquisition parameters: NS=2, DS=4, AQ=4.19.  Entry | Algorithm | $n_{\text{fev}}$ | Time taken ($\mathrm{s}$)
---|---|---|---
1 | Nelder–Mead | 29–40 | 1250–1726
2 | MDS | 34–53 | 1487–2314
3 | BOBYQA | 25–40 | 1143–1843
Table S23:  Solvent suppression 4-parameter (transmitter offset, presaturation
power, mixing time, and presaturation duration) optimisations. The POISE
routine used was: . Key acquisition parameters: NS=2, DS=4, AQ=4.19.

The best of each of these optimisations are shown in the next figure.

Figure S13:  Insets of 1D NOESY / presaturation spectra of a sample of rodent
urine, obtained before and after optimisation of up to 4 parameters. The time
required for each optimisation is indicated on each spectrum.

###  S11.4 1D NOESY pulse programme

This pulse programme is almost the same as the Bruker default noesygppr1d,
except that the presaturation power is controlled via the CNST27 parameter
(which represents the RF amplitude in $\mathrm{Hz}$).

⬇

;noesygppr1d_jy

;

;download from: https://git.io/Jnokc

;

;modified from noesygppr1d to use CNST20 as presaturation RF amplitude (in Hz)

;Jonathan Yong - University of Oxford - 20 June 2021

;

;$CLASS=HighRes

;$DIM=1D

;$TYPE=

;$SUBTYPE=

;$COMMENT=

;$RECOMMEND=

#include <Avance.incl>

#include <Grad.incl>

”d12=20u”

”plw9=plw1*(cnst20/(250000/p1))*(cnst20/(250000/p1))”

”acqt0=-p0*2/3.1416”

1 ze

2 30m

# ifdef FLAG_BLK

4u BLKGRAD

# else

4u

# endif /*FLAG_BLK*/

d12 pl9:f1

d1 cw:f1 ph29

4u do:f1

50u UNBLKGRAD

p16:gp1

d16 pl1:f1

p1 ph1

4u

p1 ph2

d12 pl9:f1

d8 cw:f1

4u do:f1

p16:gp2

d16 pl1:f1

# ifdef FLAG_BLK

4u

# else

4u BLKGRAD

# endif /*FLAG_BLK*/

p0 ph3

go=2 ph31

30m mc #0 to 2 F0(zd)

# ifdef FLAG_BLK

4u BLKGRAD

# else

4u

# endif /*FLAG_BLK*/

exit

ph1=0 2

ph2=0 0 0 0 0 0 0 0 2 2 2 2 2 2 2 2

ph3=0 0 2 2 1 1 3 3

ph29=0

ph31=0 2 2 0 1 3 3 1 2 0 0 2 3 1 1 3

;pl1 : f1 channel - power level for pulse (default)

;pl9 : f1 channel - power level for presaturation

;p0 : for any flip angle

;p1 : f1 channel - 90 degree high power pulse

;p16: homospoil/gradient pulse

;d1 : relaxation delay; 1-5 * T1

;d8 : mixing time

;d12: delay for power switching [20 usec]

;d16: delay for homospoil/gradient recovery

;ns: 8 * n, total number of scans: NS * TD0

;ds: 4

;cnst20: presaturation power (Hz)

;use gradient ratio: gp 1 : gp 2

; 50 : -10

;for z-only gradients:

;gpz1: 50%

;gpz2: -10%

;use gradient files:

;gpnam1: SMSQ10.100

;gpnam2: SMSQ10.100

;preprocessor-flags-start

;FLAG_BLK: for BLKGRAD before d1 rather than go

; option -DFLAG_BLK: (eda: ZGOPTNS)

;preprocessor-flags-end

##  S12 Diffusion NMR

All spectra in this section were acquired on a $600\text{\,}\mathrm{MHz}$
Bruker AVIII spectrometer, equipped with a Prodigy N2 broadband cryoprobe with
a nominal maximum $z$-gradient strength of
$66\text{\,}\mathrm{G}\text{\,}{\mathrm{cm}}^{-1}$. The sample used was
$45\text{\,}\mathrm{m\textsc{m}}$ andrographolide in DMSO-$d_{6}$.

###  S12.1 Sequential optimisation strategy

Traditionally, setting up a diffusion experiment is considered to be a tedious
process. At the very least, one must search for a range of gradient amplitudes
($G_{\mathrm{min}}$ through $G_{\mathrm{max}}$) which provide an appropriate
range of signal attenuation: thus, for example, the increment acquired with
$G_{\mathrm{max}}$ must exhibit a certain level of attenuation with respect to
$G_{\mathrm{min}}$. If $G_{\mathrm{max}}$ is too large or too small, then the
resulting fitting will be less accurate. In more complicated cases, the
diffusion delay $\Delta$ and/or the gradient pulse duration $\delta$ must also
be tweaked.

Here, we use the oneshot DOSY experiment (Magn. Reson. Chem. 2002, 40 (13),
S147–S152, DOI: 10.1002/mrc.1107); the process can easily be adapted to all
other DOSY variants. Throughout the optimisation, the script maintains a
“reference” 1D spectrum acquired with $G_{\mathrm{min}}=10\%$. The entire
optimisation can be run using a single wrapper Python script, dosy_opt.py.

We seek to first optimise $\Delta$ by performing stepwise increases from a
small initial value, until a separate spectrum, obtained with
$G_{\mathrm{max}}=80\%$, displays “sufficient” attenuation with respect to the
reference spectrum. Here, we define “sufficient” as 75% attenuation. To
determine this, the script leverages the dosy_aux routine and cost function.
Since this is not a traditional “optimisation” in the sense that the value can
only go upwards, we use the wrapper script to check the value of the cost
function: thus we set the maxfev option to 1, such that on every function
evaluation, control is returned to the wrapper script. The routine is (note
that lb, ub, and tol have no use in this routine and can be set to any valid
value):

⬇

{”name”: ”dosy_aux”, ”pars”: [”gpz1”], ”lb”: [79.0], ”ub”: [81.0], ”init”:
[80.0], ”tol”: [0.01], ”cf”: ”dosy_aux”, ”au”: ”poise_1d”}

The dosy_aux cost function is defined by

$f_{\text{att}}=\frac{\sum_{i}\mathbf{S}_{i}}{\sum_{i}\mathbf{R}_{i}}-0.25$
(7)

where $\mathbf{R}$ is the reference spectrum (acquired with a minimum gradient
amplitude of 10%) and $\mathbf{S}$ is the maximally attenuated spectrum
(acquired with a gradient amplitude of 80%).

Once a suitable diffusion delay has been found, we then carry out a simple
POISE optimisation to find the actual value of $G_{\mathrm{max}}$ that
provides 75% attenuation. This utilises the dosy routine and cost function.
The routine is as follows:

⬇

{”name”: ”dosy”, ”pars”: [”gpz1”], ”lb”: [20.0], ”ub”: [80.0], ”init”: [50.0],
”tol”: [2.0], ”cf”: ”dosy”, ”au”: ”poise_1d”}

and the dosy cost function is simply the absolute value of $f_{\mathrm{att}}$.

The total time taken for the entire process depends on how many times $\Delta$
must be increased from its initial value (roughly 1 additional minute per
increase). In order to minimise this, it is recommended to start with a
sensible value of $\Delta$ based on known precedent. The second step takes
approximately 3 minutes in our hands regardless of optimisation algorithm, for
an overall total time of ca. 5 minutes (using the parameters DS=2, NS=4).

In Figure S14, we plot the typical decay profiles obtained for CH and OH
protons in andrographolide, obtained from a oneshot DOSY experiment using
parameters optimised according to the procedure described above. The two
profiles are slightly different because of chemical exchange between the OH
protons and water present in the sample. This data is clearly amenable towards
further analysis using processing tools such as Bruker’s Dynamic Centre
software or the DOSY Toolbox.

Figure S14: Example diffusion decay profiles obtained using optimised values
of $G_{\mathrm{max}}=75.6\%$ (equivalent to
$49.7\text{\,}\mathrm{G}\mathrm{/}\mathrm{c}\mathrm{m}$ on the spectrometer we
used) and $\Delta=$100\text{\,}\mathrm{ms}$$ in the oneshot DOSY experiment.
$G_{\mathrm{min}}$ is $10\%$, or
$6.57\text{\,}\mathrm{G}\mathrm{/}\mathrm{c}\mathrm{m}$. Dashed lines are the
Gaussian profiles which have been fitted to the observed intensities. The
sample used was $45\text{\,}\mathrm{m\textsc{m}}$ andrographolide in
DMSO-$d_{6}$.

###  S12.2 Python script for automation

The steps to use this script are:

  1. 1.

Make sure that you have enabled the DOSY routines and cost functions. The
routines can be found in the .../py/user/poise_backend/example_routines
directory and should be copied to the .../py/user/poise_backend/routines
directory; the cost functions are in
.../py/user/poise_backend/costfunctions.py and must be uncommented before they
can be used.

  2. 2.

Download this Python script. There is a GitHub link in the source code below.

  3. 3.

Set up a 2D diffusion experiment (most easily done by reading in a parameter
set). The initial gradient amplitudes do not matter, but the diffusion delay
should be set to something that is not too large (since the script only
increases the diffusion delay if it is found to be too short).

  4. 4.

If the particular DOSY experiment has not been registered in this script yet,
then do so by modifying the dictionary near the top of the script. This only
needs to be done once. (If the DOSY experiment has never been set up before,
then you need to find (or create) a 1D version of the diffusion experiment. It
will basically be the same as the 2D version, except that the encoding
gradient is not incremented.)

  5. 5.

Run this Python script (just type dosy_opt in the TopSpin command line).

⬇

”””

dosy_opt.py

———–

Download from: https://git.io/JUHOY

Script that optimises parameters for a diffusion experiment. Currently works

with oneshot DOSY (Pelta et al., Magn. Reson. Chem. 2002, 40 (13), S147-S152,

DOI: 10.1002/mrc.1107), as well as the stegp1s and ledgp2s DOSY variants which

are standard TopSpin pulse programmes.

This is meant to be called from the TopSpin command line on a 2D diffusion

experiment.

Note that this script does not tweak the gradient lengths (”little delta”), so

these should be set to a sensible value (e.g. 1 ms) before running this
script.

The algorithm is as follows:

1. Create 1D versions of the diffusion experiment, one with the minimum

gradient (10%, the ”reference spectrum”, expno 99998) and one with the

maximum gradient (80%, ”optimisation spectrum”, expno 99999).

2. Acquires both spectra.

3. Check whether the opt. spectrum has at least 75% signal attenuation

relative to the ref. spectrum. If not, increases the diffusion delay

(”big Delta”) and repeats steps 1 and 2\.

4. If the attenuation is sufficient, it then runs a POISE optimisation on the

opt. spectrum to find the gradient strength at which 75% signal attenuation

is achieved.

5. Once this gradient strength has been found, returns to the 2D experiment

and runs the ‘dosy‘ Bruker AU programme using the new parameter.

SPDX-License-Identifier: GPL-3.0-or-later

”””

# Dictionary of DOSY pulse programmes, to be entered in the form

# ”2d_pulprog”: [”1d_pulprog”, ”GRAD_AMPLITUDE_PARAM”, ”DIFF_DELAY_PARAM”]

# In more detail:

# - 2d_pulprog is the name of the 2D DOSY pulse sequence.

# - 1d_pulprog is the name of the corresponding 1D DOSY pulse sequence, which

# should be just a single increment of the 2D sequence, but with Difframp

# removed. See the doneshot_nd_jy pulse programmes for an example.

# - GRAD_AMPLITUDE_PARAM is the parameter name in the 2D DOSY pulse sequence

# that is multiplied by Difframp. In the 1D DOSY pulse sequence it should

# not be multiplied by Difframp, so simply corresponds to the gradient

# strength. Note that you need a space between GPZ and the number.

# - DIFF_DELAY_PARAM is the parameter name in both 1D and 2D sequences that

# corresponds to the diffusion delay. Again, a space between D and the

# number is required.

# Other DOSY variants can be entered here to ”register” them and make them

# available for use with this script.

pulprog_dict = {

# doneshot_nd_jy can be downloaded from https://git.io/JUHOX and

# https://git.io/JUHOP.

”doneshot_2d_jy”: [”doneshot_1d_jy”, ”GPZ 1”, ”D 20”],

# the following two entries are Bruker standard sequences and already come

# with TopSpin.

”stegp1s”: [”stegp1s1d”, ”GPZ 6”, ”D 20”],

”ledgp2s”: [”ledgp2s1d”, ”GPZ 6”, ”D 20”],

}

def main():

# Look up the 1D pulse programme and associated parameters based on the

# current 2D pulse programme.

pp2d = GETPAR(”PULPROG”)

if pp2d in pulprog_dict:

pp1d, gpzparam, dparam = pulprog_dict[pp2d]

else:

MSG(”dosy_opt: unsupported pulse programme. New pulse programmes can”

” be added in the dosy_opt.py script.”)

EXIT()

# Create new routines. We need to do this within this script since the

# routine depends on the GPZ parameter being optimised; we can’t hardcode

# that in a predefined routine.