Motion prediction enables simulated MR-imaging of freely moving model organisms

Abstract

Magnetic resonance tomography typically applies the Fourier transform to k-space signals repeatedly acquired from a frequency encoded spatial region of interest, therefore requiring a stationary object during scanning. Any movement of the object results in phase errors in the recorded signal, leading to deformed images, phantoms, and artifacts, since the encoded information does not originate from the intended region of the object. However, if the type and magnitude of movement is known instantaneously, the scanner or the reconstruction algorithm could be adjusted to compensate for the movement, directly allowing high quality imaging with non-stationary objects. This would be an enormous boon to studies that tie cell metabolomics to spontaneous organism behaviour, eliminating the stress otherwise necessitated by restraining measures such as anesthesia or clamping. In the present theoretical study, we use a phantom of the animal model C. elegans to examine the feasibility to automatically predict its movement and position, and to evaluate the impact of movement prediction, within a sufficiently long time horizon, on image reconstruction. For this purpose, we use automated image processing to annotate body parts in freely moving C. elegans, and predict their path of movement. We further introduce an MRI simulation platform based on bright field videos of the moving worm, combined with a stack of high resolution transmission electron microscope (TEM) slice images as virtual high resolution phantoms. A phantom provides an indication of the spatial distribution of signal-generating nuclei on a particular imaging slice. We show that adjustment of the scanning to the predicted movements strongly reduces distortions in the resulting image, opening the door for implementation in a high-resolution NMR scanner.

About the Authors:

Markus Reischl

Contributed equally to this work with: Markus Reischl, Mazin Jouda

Roles Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing

Affiliation: Institute for Automation and Applied Informatics, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

Mazin Jouda

Contributed equally to this work with: Markus Reischl, Mazin Jouda

Roles Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – review & editing

Affiliation: Institute of Microstructure Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

ORCID logo http://orcid.org/0000-0002-1226-1174

Neil MacKinnon

Roles Data curation, Formal analysis, Investigation, Methodology, Project administration, Writing – review & editing

Affiliation: Institute of Microstructure Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

Erwin Fuhrer

Roles Methodology, Visualization, Writing – review & editing

Affiliation: Institute of Microstructure Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

ORCID logo http://orcid.org/0000-0002-5316-6579

Natalia Bakhtina

Roles Investigation, Software, Validation

Current address: ETH Zurich, Zurich, Switzerland

Affiliation: Institute of Microstructure Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

ORCID logo http://orcid.org/0000-0002-9856-2559

Andreas Bartschat

Roles Investigation, Software, Visualization

Affiliation: Institute for Automation and Applied Informatics, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

Ralf Mikut

Roles Conceptualization, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing

Affiliation: Institute for Automation and Applied Informatics, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

ORCID logo http://orcid.org/0000-0001-9100-5496

Jan G. Korvink

Roles Conceptualization, Data curation, Funding acquisition, Project administration, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Institute of Microstructure Technology, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

ORCID logo http://orcid.org/0000-0003-4354-7295

Introduction

A major challenge in biological science is to relate molecular regulation at the cellular level to response and behaviour at the organism level. Knowing this relationship lies at the foundation of every disease, and indeed also in understanding the healthy organism. An experiment establishing this relationship, as schematically shown in Fig 1, requires i) in vivo cellular-level detection of regulation-relevant molecules, such as metabolites and their production rates, ii) precise access to and application of biological perturbation mechanisms, and iii) an unburdened (naturally responding) organism. Clearly, even taken individually, these items are very hard to achieve in general, and correspond to major research areas in their own right.

[Figure omitted. See PDF.]

Fig 1. A freely moving organism is subject to a stimulus (mechanical, chemical, light, etc.) causing a response at all levels of detail (metabolomic, behaviour, etc.).

The experiment designer relates the correlated output to a biological hypothesis and may adapt the stimulus.

https://doi.org/10.1371/journal.pcbi.1006997.g001

Small organisms such as C. elegans are biological model organisms useful for studying many human disorders, including neurodegenerative diseases [1, 2]. Model organisms have been the mainstay of biological sciences for decades, and thus a broad knowledge base already exists starting from the genome level, through developmental cycles, and up to behavioural response to applied stress. These platforms offer the opportunity to address the connection between molecular phenotype, which can be conveniently implemented due to the very rapid yet standardized life cycle of C. elegans, to behaviour. What remains is the technological challenge of satisfying the three requirements for robustly linking phenotype to behaviour.

Nuclear magnetic resonance (NMR), which is a noninvasive and non-destructive technique, together with its imaging modality (MRI), is a strong candidate as the analytical method of choice, towards the ultimate goal of in vivo measurement of the molecular response. NMR is based on exciting the spin-active nuclei of the magnetized organism with radio frequency (RF) signals, and detecting their response via induced RF signals. Many nuclei are NMR sensitive, but in vivo molecular concentrations of metabolites are typically below millimolar levels, requiring high sensitivity to detect. NMR is a non-ionizing technique, making it superior to computed tomography (CT) based on X-rays.

Although MRI microscopy currently only achieves spatial resolutions down to 4 μm, new techniques are partially overcoming these limits, such as the use of nitrogen vacancy centers in diamond, or the use of hyperpolarisation techniques. MRI is a Fourier imaging technique that uses magnetic field gradients to selectively excite parts of the object, and perform consecutive phase and frequency encoding of the excited spins. This so-called spatial encoding makes MRI a relatively slow imaging technique, as a single MRI image requires multiple consecutive signal acquisitions from the object that is being imaged. In conventional MRI, it is assumed that the organism (or object) being imaged is fixed in space, so that spatio-temporally varying magnetic fields are only due to the technical system. Any movement of the organism results in measurements from shifted spatio-temporal positions resulting in image artifacts and ensuing difficulty in data interpretation (see Fig 2). Motion-induced artifacts are the challenge since spatially localized spectroscopy (required for molecular profiling) necessarily requires repeated measurements to bring the signal level above the noise.

[Figure omitted. See PDF.]

Fig 2. Concept of moving frame imaging: An object of interest is placed close to an NMR sensor (coil).

For conventional cartesian acquisition of an MR image, N repetitions are required to fill the k-space matrix. The signal of one repetition is used to fill one k-space line. This acquisition scheme results in a total acquisition time of TR ⋅ N. The time scales of worm motion and repetition time TR are of the same magnitude, and with a standard Eulerian or fixed frame MR procedure, strong motion blurring occurs. To enable imaging free of blurring of a freely moving worm, the concept of a Lagrangian moving frame is required.

https://doi.org/10.1371/journal.pcbi.1006997.g002

In clinical applications of MRI, free body motion remains a challenge and there is a concerted effort underway to i) collect the MRI data faster than the motion; ii) collect the MRI data at moments when motion is minimal (triggering); iii) track patient motion and correct the MRI data during post processing; iv) track the motion and guide the spatial encoding to reflect the instantaneous geometrical configuration [3–6]; v) predict the motion and adjust the spatial encoding system in real-time, e.g., the prospective motion correction based on respiratory motion prediction [7]. Whilst these methods have been highly successful to control artefacts due to breathing, heartbeat, and low amplitude head movement, their assumptions for the kinematics of the underlying movement or the periodicity of motion are limiting. For example, the head is motion-captured in a model that assumes rigid six-degree-of-freedom body movement involving translations and rotations (x, y, z, θx, θy, θz) along the three orthogonal Cartesian axes. The direct translation of these techniques to MR measurements of small samples is not straightforward primarily because of the reduced size, more complex organism motion including writhing and wiggling, and rapid displacement across the detector’s sensitive region. Methods used to immobilise an organism can be considered to avoid these motion artifacts, such as clamping or freezing; however, such drastic measures typically introduce an undesired stress response into the molecular profile. The challenge of molecular measurement of non-stressed, small model organisms therefore still remains open. Given the advances in image processing, we believe there is an opportunity to address this challenge computationally.

Advances in computer vision have revolutionized the speed and accuracy with which image analyses can be performed, and aim to reduce the need for expert knowledge e.g. in medical image analysis [8] and automated experiment handling [9]. These technologies are entering public awareness through the automation of highly complex processes, such as trajectory generation for self-driving vehicles (road, aerial) and surveillance for public safety. This is accomplished by real-time processing of dynamic images, in which accuracy and speed are of paramount importance. To name an example, advances in image and data processing algorithms are expected to make real-time dynamic spectrum imaging (achieving a hyperspectral imaging cube or hypercube) possible at all electromagnetic wavelengths. For instance, recording image and spectral data over 500 × 500 pixels, with a spectral resolution of 5nm over the visible spectrum, at 5 f s−1 has already been demonstrated [10]. Enormous data storage and processing power is required to perform such operations with sufficient spectral, spatial, and temporal resolution, and has motivated an increased effort in sparse imaging modalities.

MRI also suffers from limitations in acquisition speed, which in this case are due to physical constraints such as the relaxation times. This renders MRI a relatively slow imaging technique compared to techniques that are based on the visual range of wavelengths. In MRI, data is acquired in -space ( being the spatial frequency), and the observer requires the transformed data in space. Hence, a large number of acquired data is associated with each final image pixel or voxel, and motion or geometrical warping during imaging will introduce errors into the reconstruction algorithms. To obtain sufficient acquired signal power per voxel volume, object tracking must be implemented in order to ensure correct spatial and spectral co-localization over time. This is especially important in MRI microscopy, which depends on accumulated sampling for sufficient image resolution. With the novel capabilities of image and data processing algorithms, real-time image processing becomes feasible. This enables the extraction of information and motion prediction based on a conventional video stream which then can be used to steer the MR sequence in real-time. The prediction horizon will therefore be dependent on the signal acquisition time at a given frame rate and voxel spatial resolution.

In this contribution we consider the preconditions for performing MRI experiments on unburdened small organisms. We will focus our attention on the nematode C. elegans, mainly because it fits in best with our own efforts towards in vivo metabolomic profiling, but we will address the underlying problem in more generality so that it is relevant also for other organisms.

In this contribution we explore the possibility to completely remove the requirement of organism immobilisation, by providing the host observing technical system (e.g., a spectrometer or microscope) with a real time co-evolving Lagrangian coordinate system centered in the organism that provides the organism’s current center-of-gravity position and shape, paired with a robust prediction of the organism’s future position and deformed shape.

There is a multitude of publications of so-called motion-trackers, which mainly quantify rigid body motions and ignore strain fields. Important ones are [11] which categorises behavior and morphology features out of C. elegans-videos by applying segmentation and tracking algorithms and [12], which applies methods of machine vision, data processing and tracking to evaluate drug assays and [13] to follow multiple worm within the same environment. Furthermore, motion decomposition methods using so-called Eigenworms apply methods from oscillation analysis to divide the worm motion into so called eigenmodes [14]. There are automated software-packages being able to quantitatively assess movement parameters (e.g. Track-A-Worm [15], Parallel Worm Tracker [16] etc.). However, there is a lack of an easy-to-use and sufficiently fast prediction algorithm for individually chosen positions within a single worm to be used to adjust the gradient system following a body part to be imaged. The prediction horizon available for the host system depends on the speed of image processing and latency of the capturing hardware (essentially the frame rate), as well as the smoothness of the organism’s motion. To gain insight into how movement prediction can enhance MRI signal detection, the paper introduces:

* Data of a virtual phantom, combining high resolution electron microscope slice images [17] with conventional video-recordings of C. elegans moving in a Petri dish (Section Phantom Generation);

* A new concept for location prediction in C. elegans, and detailed movement prediction, exploiting characteristic worm movements (Section MRI Simulation);

* A computational platform, being able to evaluate the outcome of MR imaging with and without adapting the imaging parameters based on the movement prediction (Section MRI Simulation).

* A measure for simulation success that compares prediction to the ‘true’ image through a similarity measure sxy, a technique first introduced in [18] and described in detail in the supporting information.

Using this simulation model, we demonstrate the capabilities of the motion-prediction algorithm in MRI.

Materials and methods

Phantom generation

The study is based on eight 10 s duration AVI-videos (sampling frequency 12 Hz) of C. elegans. The recording was done in a controlled environment with a fixed camera and constant illumination over time. The worm was enclosed in a technical setup including a microfluidic channel and a Petri dish of the host system, suitable for optical recording and other real-time measurements. Typically, if natural state studies are relevant, the organism can also be provided with an optically transparent nutritional substrate such as a gel containing E. coli bacteria.

To simulate MR imaging of the moving worm, we artificially linked transmission electron microscope (TEM) images to the video of the worm, such that virtually scanning each location within the worm would deliver a simulated high-resolution MR-image. At each time-frame, the 1.2mm worm was segmented into 50 slices perpendicular to its center line. This resulted in a reasonable slice thickness of roughly 24μm. Subsequently, we assigned voxel MR signals adapted from the TEM images to each slice as follows: first, we removed the background from the TEM images and scaled them such that the dimensions of the body part in each slice were realistic. Second, we inverted the color map of the images, since the white regions of the TEM images corresponded to low density material, and therefore would appear dark in MR imaging due to the low proton signal intensity. This color inversion is just to make the sample images look like MR images. The reader should bear in mind that TEM and MRI use completely different imaging mechanisms, and therefore converting one to the other is not a straightforward task. In fact, having true MRI sample images is not essential to prove the efficiency of motion prediction in imaging enhancement, and indeed any arbitrary set of morphologically correct sample images would work. Finally, we reduced the number of pixels of the TEM images to 64 × 64, which corresponded to an MRI in-plane resolution of approximately 1.6 μm, and then assigned the voxel signals to the virtual MR slices. Although such a high volumetric resolution lies beyond the capabilities of currently available MRI scanners, it was chosen on purpose to allow more accurate assessment of the prediction algorithm. It is thus a good compromise between being close to a realistically achievable resolution, and being high enough for the structural similarity measure to work more accurately. Furthermore, considering the rapid advances in the field of MRI including higher fields, stronger gradients, and hyperpolarization techniques, the mentioned resolution of 1.6μm is believed to be achievable in the near future.

Motion prediction

New concept.

In order to adjust the gradient system of the MRI, the worm needs to be detected and its future positions need to be predicted. Therefore, we introduced a new concept of real-time image processing of the worm, suitable for any video format, regardless of color model, and also independent of background structures.

Fig 3a shows a block diagram of the basic steps of the concept. Starting with a raw image matrix Xraw, a robust estimation of the worm and its position Xworm within the video were determined. Based on the segmented worm, characteristic points such as center of gravity (COG) and position of the head region were determined. To introduce a coordinate system along the center line of the worm, a skeletonization was applied, see Fig 3b. Assuming that body regions excluding the head would move along the center line, the worm velocity was used to predict positions xc,p of arbitrary selected points xc within the worm. A detailed description is given in the supporting information.

[Figure omitted. See PDF.]

Fig 3.

a: Concept for the prediction of defined locations within a worm, b: Parametrization of the worm, c-j: Snapshot of a sample video. c: grayscale image/video Xraw[k], d: Image with subtracted background, e: difference image, f: overlay with worm, g: skeletonization, h-j: Point of interest (POI, gray marker) and predictions for Δk = 30 (= 2.5 s) time samples (white marker). h-j: Time samples k = 50, k = 80, k = 110. Figures show the position every 30 samples and the prediction of the POI 30 samples ahead. Thus, the similarity of POI(80) vs. POIpred(80) and POI(110) vs. POIpred(110) shows the quality of the prediction.

https://doi.org/10.1371/journal.pcbi.1006997.g003

Worm detection.

This section briefly covers the worm detection and movement prediction. The preprocessing for the worm detection aims to reduce the computational complexity and consists of a grayscale conversion (standard Matlab conversion) of the videos as well as a decrease of the resolution by a factor of 9 (Fig 3c) (Decreasing the resolution proves real time processing to be possible / usage of low-end imaging hardware to be accurate (see Discussion)). For the segmentation of the worm, the background of the video is estimated and removed from the images (Fig 3d).

The direction of motion is determined by relying on the position of the head and the COG. The COG is calculated based on the foreground pixels of the segmentation and is smoothed to obtain variations for future predictions (1st order low-pass filtering over time). The position of the head is computed utilizing difference-images (Fig 3e) of the last 10 frames in the video. Fig 3f shows the detected head.

Coordinate prediction.

The position of the COG moves linearly and is predicted by linear extrapolation based on its past five positions. To predict the position of an arbitrary point of interest (POI) within the worm (which can easily be identified e.g. by a mouse click in the an NMR-simulation), we use the skeleton line as the center line of the segmentation and introduce a normalized coordinate s along it (head: s = 0, tail: s = 1, see Fig 3b and 3g).

Assuming that each worm segment moves according to the current shape of the worm, following its predecessor segment with the velocity of the worm (The assumption is not perfectly valid, but it significantly simplifies processing and delivers reasonable results (see also Fig 4)), the velocity and the shape is used to predict the location sc of the POI after Δk time steps using sc = s − vΔk, (Fig 3h–3j) (If the prediction horizon is chosen too high, negative sc are avoided by setting the prediction to the topmost point (s = 0)).

[Figure omitted. See PDF.]

Fig 4. Evaluation (mean values and standard deviation) of the dataset given in section phantom generation.

a: fixed location in the worm (s = 0.9), varying prediction horizons, b: fixed prediction horizon (Δk = 12, 1 second, varying locations, prediction horizon too high for s < 0.3).

https://doi.org/10.1371/journal.pcbi.1006997.g004

The prediction quality is evaluated using an Euclidean distance. Given the prediction horizon Δk, the predicted position based at time point k + Δk is calculated (x1,c,p[k + Δk], x2,c,p[k + Δk]) and the distance to the true position (x1,c[k + Δk], x2,c[k + Δk]) is measured. The average for all time samples is termed the mean prediction error. The dependency of the prediction accuracy to the prediction horizon, as well as to the predicted position of the worm, is given in Fig 4: The prediction error increases if the prediction horizon becomes larger. Regions near the head cannot be reliably predicted, since the head moves at a much higher frequency and in random patterns while scanning the surroundings and deciding on the moving direction. Furthermore, if the prediction horizon is set too high the predicted point lies outside the current shape of the worm (e.g. for Δk = 12 for all s < 0.3) (depending on the speed of the worm, however wild type worms mostly have the same velocity of approximately 0.13 mm s−1) and the topmost point (s = 0) is chosen as prediction, resulting in deviations).

For all other regions, the uncertainty of the prediction stays roughly constant.

Results

MRI simulation

Simulation-platform.

The magnetic resonance imaging process is simulated using Matlab. The worm is assumed to reside within a strong and constant magnetic field B = (0, 0, Bz) of the MRI scanner. The software mimics a simplified yet acceptably accurate MR image acquisition process via a standard gradient-echo sequence [19] that takes place during worm motion as follows:

* An image slice at a parametrised position s along the worm axis is defined as the image plane perpendicular to the worm axis. As previously mentioned, the worm is segmented into 50 slices whose voxels’ intensities are obtained from the TEM images. Thus s is the center of one of those slices. Ideally, to do MRI of a slice s, the center of the gradient system (the field of view FoV), GC, should be exactly at s over the entire signal acquisition process. This is unfortunately not the case as the prediction imposes some error resulting in GC being at another slice, at the desired slice but off-center, or both. In this case the MR excitation is done by calculating the signal intensity of each voxel in the FoV depending on the position of GC. If, for instance, GC is off-center, then some voxels of the true slice will be out of the FoV and the corresponding ones in the FoV will then be filled with zeros. The excitation step is followed by a phase encoding step, in which each voxel of the FoV is given a phase proportional to its position along one cross-sectional axis ξ. After a time delay TE/2, where TE (echo time) is the time from excitation to the center of the MR signal (echo), a frequency encoding gradient is applied, whereby each voxel of the FoV is assigned a frequency proportional to its position along an axis η perpendicular to the phase encoding axis (so that ξ ⋅ η = 0). The superposition of all voxel RF signals are simultaneously detected by a virtual coil. This so-called echo is recorded at instant TE, resulting in one line of k-space. The detection system is assumed to have a uniform spatial sensitivity.

* After a time delay TR (repetition time), which is the time it takes to repeat the sequence in order to acquire a new line of the k-space, the algorithm calculates the anticipated new position and orientation of the selected slice (as the worm would have already moved to a new position). Now the slice at this new position is excited, phase encoded, and frequency encoded, resulting in a new line of k-space being filled.

* The entire process is repeated until all the lines of k-space (in our case 64 lines) are complete.

Once the imaging procedure is complete, the program reconstructs the anticipated MR image from the k-space data via Fourier transform.

Algorithm 1 Gradient-echo MR imaging

1: procedure GRE(TE, TR)      ⊳ Echo time and repetition time.

2:  s = 0      ⊳ s parametrises slice position along worm axis.

3:  while s ≠ 1 do           ⊳ cycle through the slices.

4:   s ← s + δs                 ⊳ Next slice.

5:   t = 0                ⊳ Initialize time counter.

6:   N = 64               ⊳ Number of k-space lines.

7:   k = −π                 ⊳ Set initial phase.

8:   J = 1                     ⊳ Set index.

9:   while t ≠ N ⋅ TR do

10:    Pos. Pred.            ⊳ Predict next position.

11:    Mov. Grad.         ⊳ Update imaging gradient.

12:    Exc.                    ⊳ Excite slice.

13:    t ← t + TE/2      ⊳ Advance time to half echo time.

14:    Ph. Enc. (k + J ⋅ 2π/N)            ⊳ Phase encode.

15:    t ← t + TE/2        ⊳ Advance time to echo time.

16:    Freq. Enc.             ⊳ Frequency encode.

17:    Echo                ⊳ Signal acquisition.

18:    t ← t + TR − TE        ⊳ Advance time to next line.

19:    J ← J + 1

20:   end while

21:  end while

22: end procedure

The simulation is integrated into an interactive graphical user interface (GUI) to enable the user to set and change the imaging parameters easily and execute the processing without any programming skills (Fig 5).

[Figure omitted. See PDF.]

Fig 5. Graphical user interface of a Matlab program that simulates the MR imaging experiment of a moving worm.

The microscope movie of the worm is shown to the bottom left, and emulates real-time observation. The slice position is also shown. The top right shows the phantom MR image of the worm at the position of the slice.

https://doi.org/10.1371/journal.pcbi.1006997.g005

The simulated MRI signals can be chosen by the user and are either MR-like signals translated from real transmission electron microscopy (TEM) images, or any virtual MR-like images provided by the user.

To start the simulation, an arbitrary slice of interest (with coordinate s) within the worm’s body is selected. For this point, the simulated MR image of a slice without movement is shown. In the display, the slice position s is denoted by a red line perpendicular to the center line of the worm, while the predicted position to which the virtual imaging gradient coordinate is set is denoted by a green line.

Furthermore, the software allows the user to set the TR and TE imaging parameters, and to choose the prediction horizon (a value between 1 and 10 frames), which refers to the number of frames ahead for which the predicted worm position will be calculated. Our algorithm does not correct for the worm motion occurring during recording of a line of k-space. This is equivalent to the assumption that the MRI pulse sequence is based on very short echo times, which is possible, but taken at the expense of increased acquisition bandwidth and thus reduced signal-to-noise ratio (SNR).

Simulation paradigms.

We performed three MRI simulations using the optical video dataset from Section Phantom Generation to:

1. confirm that the prediction enhances the MR imaging of the moving worm in general (Simulation 1);

2. measure the efficiency of the prediction algorithm as the desired resolution of the MR images increases in comparison to the resolution of the optical video used for the prediction (Simulation 2); and

3. measure the effect of an increased prediction horizon (the number of frames ahead for which the position is predicted, Simulation 3).

In all the simulations, the repetition time (TR) was set to 83 ms, which corresponds to the frame rate of the optical video on which the prediction was based.

We evaluate the simulation by quantifying the structural similarity between the true image and simulated image as introduced in [18]. Identical images return sxy = 1, whereas structural inequality delivers sxy = 0.

Body position (Simulation 1).

In Simulation 1, three slices (slice definition from Section Phantom Generation, first slice from the head (s = 0.02), second from the middle body (s = 0.5), third from the tail (s = 1) of the worm) are selected and an MR imaging simulation is performed (prediction horizon Δk = 1 (83 ms)).

The results are illustrated in Fig 6a, each row shows the true slice on the left, the simulated MR image based on the prediction algorithm in the middle, and the simulated image when no position prediction is involved on the right.

[Figure omitted. See PDF.]

Fig 6.

a: effect of position prediction on the quality of the MR images in the head (s = 0.02), body (s = 0.5) and tail (s = 1) regions. The prediction horizon in this case is 83 msec. b: dependence of the MR imaging quality (structural similarity) on the position of the slice. The figure displays the simulation results from eight videos of different worms.

https://doi.org/10.1371/journal.pcbi.1006997.g006

Clearly, the prediction algorithm significantly reduces the motion artifacts that would otherwise occur if the gradient system did not follow the worm as it moves. Moreover, we observe that prediction quality varies along the worm’s body. More specifically, the prediction performs better for the slices from the body and tail when compared with prediction of the head. This is axiomatic, since the worm rapidly moves its head laterally whilst scavenging for nutrition, thus the motional entropy of the head is higher (see also Fig 4). Fig 6b shows the simulation results of eight videos of different worms. The abscissa represents the slice position starting from the head (s = 0) to the tail (s = 1), while the ordinate shows the structural similarity between the prediction-based MR images and the original slices. To very good agreement with Figs 4 and 6b shows that the quality of the prediction-based MR images are higher for slices from the body and tail than the images taken near the head where the prediction uncertainty is usually higher.

Image resolution (Simulation 2).

In Simulation 2, a slice from the middle body section of the worm (s = 0.72) was selected and an MRI simulation was performed for different resolution ratios (= MR image: optical video). The idea here is to explore the effect of changing the MRI resolution (voxel size) on the quality of the prediction-based imaging. Imagine, for example, that one pixel of the optical video is 10 μm × 10 μm and the desired MRI isotropic resolution is also 10 μm × μm × 10 μm × 10 μm, then one pixel error in prediction results in one voxel error in the MRI k-space. On the other hand, if higher MRI resolution is wished, then one pixel error in prediction results in more than one voxel error in the MRI k-space. Decreasing the optical resolution can be used to speed up image processing, and thus decreasing the prediction horizon, if needed. The ratios of 4:1 and 2:1 are heuristically chosen, for which each pixel in the optical video respectively corresponds to 4 and 2 pixels in the MR image. Fig 7 demonstrates the results of this simulation—the rows correspond to the different resolution ratios while the columns depict the original slice, the image with prediction, and the image without prediction. The figure shows that the efficiency of the prediction declines as the desired resolution of the MR image increases (pixel size decreases), or alternatively, the optical resolution should be higher or at least equal to the desired MR resolution for a high quality prediction-based MR image. Indeed, the commercial optical imaging solutions can easily meet such demands of high resolution; however, upon implementation, any increase in number of pixels will be at the expense of the prolonged processing time needed for prediction. Nevertheless, this can be overcome by utilizing more powerful processors and by employing techniques such as parallel computing.

[Figure omitted. See PDF.]

Fig 7. Effect of decreasing the ratio of the MR image resolution to the optical video resolution on the imaging quality for a prediction horizon of 1.

The rows show the results for MR image:video resolution ratios of 4:1 and 2:1 respectively. Each row displays, from left to right, the reference slice, the MR image with position prediction, and the MR image without position prediction.

https://doi.org/10.1371/journal.pcbi.1006997.g007

Prediction horizon (Simulation 3).

Simulation 3 varies prediction horizons for one slice (s = 0.52). Fig 8 illustrates the results of this simulation: Fig. (1a–10a) are the prediction-based simulated images for prediction horizons of Δk 1 to 10, respectively. In contrast, Fig. (1b–10b) display the results when prediction is not in action. Because only a few lines of the k-space are collected from the correct slice (depending on TR and the speed of the worm), the images in this figure exhibit a noticeable decrease in quality as the prediction horizon increases, leading to the conclusion that one should, whenever possible, minimize the prediction horizon. Of course, in an actual hardware realization of the system, the choice of the prediction horizon will be bounded by the speed of the image acquisition and processing units.

[Figure omitted. See PDF.]

Fig 8. Effect of enlarging the prediction horizon on the quality of the imaging (s = 0.52, see Fig 6).

The images (1a-10a) show the reconstructed MR images based on the prediction algorithm for prediction horizons from 1 to 10 frames. The images (1b-10b) show the reconstructed MR images when no prediction is involved and when the initial position of the imaging gradient is set to the position of the slice after 1 (83 ms) to 10 (830 ms) frames respectively. (c) The similarity measure of the simulated image versus the prediction horizon.

https://doi.org/10.1371/journal.pcbi.1006997.g008

The effect of increasing the prediction horizon is described quantitatively in Fig 8c by measuring the structural similarity between the original slice and the simulated image for both cases with and without prediction. The results confirm the variation in predictability along the length of the worm. The prediction horizon degradation is linear for the rear two-thirds of the worm, but falls off more rapidly for the head section, as expected from Fig 4.

Moreover, a statistical assessment of the effect of the prediction horizon on the quality of the MR imaging was done: Regarding simulation results of four slices along the worm using the given eight videos it can be shown that the accuracy of prediction and thus the MRI quality decays with the increased prediction horizon but is in mean for all parameter combinations roughly three times better than without prediction.

Discussion

The results demonstrate that our prediction algorithm can markedly improve MR image quality of arbitrarily moving and deforming objects. The complete processing pipeline (including functions which will not be used in real-time processing, e.g. rotation, prediction of COG, bounding boxes etc.) takes 76 msec per frame on an average Laptop PC (core I7, fifth generation, 16 GB RAM) and has various possibilities for optimization (modern hardware, software parallelisation, hardware-filtering etc.). Decreasing the resolution of the optical video offers an additional option to reduce processing time up to a factor of 8 without significant loss of quality. Based on the current implementation, real-time processing is already possible for sample frequencies smaller than 13 Hz and application of the algorithm to a real MRI system is feasible.

For the predictive information to be useful, the characteristic time (velocity) of sample motion tmotion and predictive accuracy needs to be determined. This can be done by recording a freely moving sample with the desired optical setup followed by an image analysis routine as described in this report. Models of the motion kinematics can be tested in order to maximize the prediction horizon given a user-specified predictive accuracy threshold. The threshold can be selected based on the expected error introduced within a single voxel with pre-selected dimensions. The predictive horizon together with the calculation time tpred then can be used to determine if the sample is ‘MR imageable’, i.e. by comparing the timescales to those required for MRI. For example, in this report it was observed that a prediction horizon of 83 ms requiring a calculation time of 76 ms yielded a predictive accuracy of approximately 57% as indicated by the structural similarity measure of the prediction-based MRI simulation, Fig 8c. On the other hand, a prediction horizon of 830 ms requiring a calculation time of 76 ms resulted in a prediction accuracy of approximately 30%.

With the timescales of motion and predictive calculation defined, one must evaluate whether MRI is possible by comparing to the instrumental timescale (Only a standard cartesian MRI sequence is considered. The authors acknowledge the existence of more advanced sampling schemes, but these are outside the scope of the present discussion). The shortest relevant timescale is the time between spatial encoding steps, which in the case exemplified here is the repetition time TR (TR includes TE and the data acquisition time). It is during TR that prediction and hardware adjustments must be done prior to the subsequent spatial encoding step. In the MRI simulations described here, TR was 83 ms (including a TE of 4 ms) while tpred was 76 ms. Given the organism motion and hardware/experiment timescale regimes, one can estimate the potential for sample imaging with correction, summarized as follows:

* tmotion < TR: the object is not MR imageable without motion artifacts. Conditions to slow the natural motion of the sample should be identified and implemented.

* TR < tmotion, tpred: the object is MR imageable. Careful choice of TE and TR must be done so that the prediction calculation is complete before the next spatial encoding period. This places a restriction on the types of contrast that can be implemented.

* tpred < TR < tmotion: the object is MR imageable. There is no restriction on the contrast weightings that can be implemented.

To further improve the predictive quality, it is important to have kinematic models appropriate to the organism to be imaged. C. elegans is a convenient model for this reason given the advanced studies about behavioural phenotypes [11, 12, 20] and motion decomposition using so-called Eigenworms [14]. Extension of these models to organisms featuring similar motion characteristics should be straightforward (i.e. oscillation/undulation along the long axis—worms, snakes, swimming fish). As the kinematics becomes more complicated and/or sporadic, it will become necessary to introduce a method to guide the organism in order to introduce a predictive nature to its motion (food source, temperature gradient, etc.). Alternatively, shorter prediction horizons can be targeted together with faster MR imaging sequences such as the echo planar imaging (EPI), most likely at the expense of spatial resolution.

To bring the presented method of in-situ real-time video capturing and motion prediction experiments for MR microscopy into daily routine, we so far lack available MR-compatible optical microscopes for high-field vertical bore magnets. However, there is strong progress to reach this technical integration in the near future. Initial results have already been presented in the early 2010s and a patent describing the adaptive objective to be integrated into an optical system has been granted in 2019 under the title “MR-compatible microscope, EP 2824471B1 (see e.g. https://worldwide.espacenet.com)”. Once the remaining hardware challenges are overcome the presented procedure for a tracking control loop will give the possibility to fully exploit and expedite these novel techniques for conducting MR microscopy of free-moving, undisturbed microscopic organisms.

Supporting information

S1 Text. This file describes the principal structure of the motion prediction algorithm.

https://doi.org/10.1371/journal.pcbi.1006997.s001

(PDF)

S1 Fig. Structural similarity versus prediction horizon for slices (s = 0.1, s = 0.3, s = 0.5, and s = 0.9) for eight videos.

https://doi.org/10.1371/journal.pcbi.1006997.s002

(TIF)

Acknowledgments

We sincerely thank Dr. David H. Hall of Albert Einstein College of Medicine in New York for providing the high resolution images of C. elegans that we used as phantoms for high resolution NMR images. These lent our images more realism. The original images are from the MRC/LMB C. elegans Archive now curated by the Hall lab. These were generously donated by John White and Jonathan Hodgkin to the Hall lab.

Citation: Reischl M, Jouda M, MacKinnon N, Fuhrer E, Bakhtina N, Bartschat A, et al. (2019) Motion prediction enables simulated MR-imaging of freely moving model organisms. PLoS Comput Biol 15(12): e1006997. https://doi.org/10.1371/journal.pcbi.1006997

References

1. Bakhtina N, Korvink JG. Microfluidic Laboratories for C. elegans Enhance Fundamental Studies in Biology. RSC Advances. 2014;4(9):4691–4709.

2. Bakhtina N, MacKinnon N, Korvink JG. Advanced Microfluidic Assays for C. elegans. Online Publication—ResearchGate. 2017;.

3. Herbst M, Maclaren J, Weigel M, Korvink JG, Hennig J, Zaitsev M. Prospective motion correction with continuous gradient updates in diffusion weighted imaging. Magnetic Resonance in Medicine. 2012;67(2):326–338. pmid:22161984

4. Maclaren J, Herbst M, Speck O, Zaitsev M. Prospective motion correction in brain imaging: a review. Magnetic Resonance in Medicine. 2013;69(3):621–636. pmid:22570274

5. Herbst M, Maclaren J, Lovell-Smith C, Sostheim R, Egger K, Harloff A, et al. Reproduction of motion artifacts for performance analysis of prospective motion correction in MRI. Magnetic Resonance in Medicine. 2014;71(1):182–190. pmid:23440737

6. Godenschweger F, Kägebein U, Stucht D, Yarach U, Sciarra A, Yakupov R, et al. Motion correction in MRI of the brain. Physics in Medicine and Biology. 2016;61(5):R32–R56. pmid:26864183

7. Song H, Ruan D, Liu W, Stenger VA, Pohmann R, Fernández-Seara MA, et al. Respiratory motion prediction and prospective correction for free-breathing arterial spin-labeled perfusion MRI of the kidneys. Medical Physics. 2017;44(3):962–973. pmid:28074528

8. Li Z, Zhang X, Müller H, Zhang S. Large-scale retrieval for medical image analytics: A comprehensive review. Medical Image Analysis. 2018;43:66. pmid:29031831

9. S Y, J NB. Autonomous Injection of Biological Cells Using Visual Servoing. Experimental Robotics VII. 2001;271:169–178.

10. Ma C, Cao X, Wu R, Dai Q. Content-adaptive high-resolution hyperspectral video acquisition with a hybrid camera system. Optics Letters. 2014;39(4):937–940. pmid:24562246

11. Yemini E, Jucikas T, Grundy LJ, Brown AEX, Schafer WR. A database of C. elegans behavioral phenotypes. Nature Methods. 2013;10(9):877–879. pmid:23852451

12. Buckingham SD, Partridge FA, Sattelle DB. Automated, high-throughput, motility analysis in Caenorhabditis elegans and parasitic nematodes: Applications in the search for new anthelmintics. International Journal for Parasitology: Drugs and Drug Resistance. 2014;4(3):226–232. pmid:25516833

13. Fontaine E, Burdick J, Barr A. Automated Tracking of Multiple C. Elegans. In: 2006 International Conference of the IEEE Engineering in Medicine and Biology Society; 2006. p. 3716–3719.

14. Stephens G, Johnson-Kerner B, Bialek W, Ryu W. Dimensionality and Dynamics in the Behavior of C. elegans. PLoS Computational Biology. 2008;4(4):e1000028. pmid:18389066

15. Wang SJ, Zhao WW. Track-A-Worm, An Open-Source System for Quantitative Assessment of C. Elegans Locomotory and Bending Behavior. PLoS ONE. 2013;8.7:e69653. pmid:23922769

16. Ramot D, Johnson BE, Berry TL, Carnell L, Goodman MB. The Parallel Worm Tracker: A Platform for Measuring Average Speed and Drug-Induced Paralysis in Nematodes. PLoS One. 2008;3(5):e2208. pmid:18493300

17. Altun ZF, Herndon LA, Wolkow CA, Crocker C, Lints R, Hall DH. WormAtlas; 2002-2019.

18. Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. Image quality assessment: from error visibility to structural similarity. IEEE Transactions on Image Processing. 2004;13(4):600–612. pmid:15376593

19. Fong W. Handbook of MRI pulse sequences. Medical Physics. 2005;32(5):1452–1452.

20. Roussel N, Morton CA, Finger FP, Roysam B. A computational model for C. elegans locomotory behavior: application to multiworm tracking. IEEE Transactions on Biomedical Engineering. 2007;54(10):1786–1797. pmid:17926677