ChronoMID - Cross-Modal Neural Networks for 3-D Temporal Medical Imaging Data

TL;DR

This paper introduces ChronoMID, a cross-modal neural network approach for analyzing 3-D temporal medical imaging data, achieving high accuracy in classifying bone disease in mice by leveraging temporal features.

Contribution

It applies cross-modal CNNs to medical imaging, incorporating temporal information via timestamps and difference images, demonstrating significant accuracy improvements.

Findings

Models with temporal features outperform baseline CNNs.

Top model achieves 99.54% accuracy on unseen data.

Training on over 75,000 images validates effectiveness.

Abstract

ChronoMID builds on the success of cross-modal convolutional neural networks (X-CNNs), making the novel application of the technique to medical imaging data. Specifically, this paper presents and compares alternative approaches - timestamps and difference images - to incorporate temporal information for the classification of bone disease in mice, applied to micro-CT scans of mouse tibiae. Whilst much previous work on diseases and disease classification has been based on mathematical models incorporating domain expertise and the explicit encoding of assumptions, the approaches given here utilise the growing availability of computing resources to analyse large datasets and uncover subtle patterns in both space and time. After training on a balanced set of over 75000 images, all models incorporating temporal features outperformed a state-of-the-art CNN baseline on an unseen, balanced validation set comprising over 20000 images. The top-performing model achieved 99.54% accuracy, compared to 73.02% for the CNN baseline.

Tables

Table 1. Table 1: Accuracy scores and per-epoch training times for the different neural network classifiers, after 5 epochs and 10 epochs of training, together with the number of parameters in each model.

Method	Accuracy / %		Mean Time / s		Parameters
	5 Epochs	10 Epochs	5 Epochs	10 Epochs
CNN	$73.02$	$66.08$	$9677$	$9676$	$863883$
CNN, timestamps	$79.54$	$68.19$	$9631$	$9621$	$863947$
X-CNN, abs. diff.	$97.99$	$95.07$	$42033$	$44767$	$1643323$
X-CNN, rel. diff.	$84.28$	$85.36$	$42334$	$46816$	$1643323$
X-CNN, timestamps & abs. diff.	$98.39$	$99.54$	$44196$	$43560$	$1643387$
X-CNN, timestamps & rel. diff.	$83.43$	$75.15$	$43990$	$45356$	$1643387$

Table 2. Table 2: Percentages of observations in the “wild type” and PTH treatment classes and total number of observations, for the training, validation, and test sets, both when using image differencing and when using original images only.

		“Wild” / %	PTH / %	Total
No diff.	Train.	49.35	50.65	76663
	Val.	49.70	50.30	8519
	Test	50.24	49.76	21932
Diff.	Train.	49.46	50.54	75004
	Val.	49.34	50.66	8334
	Test	50.04	49.96	21400

Equations

$$D_{ij}=I_{ij}^{(t)}-I_{ij}^{(t_0)}$$

Full text

1 Introduction

Temporal information is especially important in the medical domain, wherein diseases develop over extended periods, treatment efficacy cannot be evaluated at the moment of administration, and temporal proximity between diseases suffered by a patient affects treatment suitability and patient recovery. In particular, studying bone remodelling (i.e. the changes over time of bone properties regulated by the activity of the bone cells) is fundamental for better understanding the effect of musculoskeletal diseases, such as osteoporosis and osteoarthritis, and to optimise the related treatments. Animal models, in particular mice, are used for studying the effect of new treatments on bone remodelling using high-resolution medical images and micro computed tomography (µCT). From µCT images, proper assessment of the morphometric parameters of the bone microstructure can be used to measure the effect of diseases or interventions [2]. Furthermore, 3D µCT images collected in vivo can be converted into biomechanical computational models (i.e. finite element models based on partial differential equations) for the non-invasive assessment of the bone’s mechanical properties [10]. Medical data, medical imaging data especially, contains a wealth of information which could, for example, inform diagnoses or detect early-stage diseases, yet remains largely untapped by automated processes. Due to the sheer quantities of data involved in many medical problems and the potential subtlety of informative patterns within it, the task of uncovering and extracting such patterns often falls outside the scope of human capabilities. Data-driven, automated approaches are essential to deriving usable insights from large medical datasets, such as those generated by scanners or ECGs. Prior work to incorporate temporal information for healthcare has focused on mining electronic patient health records for elucidating textual content or data in structured fields, thereby neglecting the temporal aspects of imaging data [3]. In this paper, we build on the successes of cross-modal convolutional neural networks [16] to incorporate both spatial and temporal information explicitly into models. We present five models, each using an alternate description of temporal information, and compare their performance against a state-of-the-art conventional convolutional neural network (CNN) for the task of bone disease classification in mice. Our approaches involved relative timestamps, corresponding to the week the scan was taken since the study began, absolute and relative difference images, and the combination of timestamps and difference images. Two groups of mice were studied: a control group of healthy mice and a group that were treated with parathyroid hormone (PTH), which can have anabolic effects leading to accelerated bone formation when administered correctly. The data for these mice — treatment status and µCT scans of the right tibia of each mouse — were provided by Dr. Enrico Dall’Ara’s team at the University of Sheffield’s Department of Oncology and Metabolism. Whilst different variants of recurrent neural networks (RNNs) are increasingly common for tasks involving sequences of data, they have limitations and drawbacks we seek to overcome. Foremost, RNNs are designed to process sequential inputs, whereas medical data may have missing data points, e.g. due to patients missing appointments, or data points provided out of order, e.g. due to different administrative bodies holding patient records because people move locations or switch healthcare providers. It is preferable to make diagnoses as data becomes available rather than waiting until it has been provided in its entirety. Furthermore, RNNs are computationally expensive, partially due to the multiple gates and units involved in memory cells, and partially due to intra-layer or backwards as well as forwards connections, leading to potentially very many parameters, in turn increasing the amount of data required to train, thus training times. Moreover, by decoupling the spatial and temporal information of the image set, we allow the networks to determine the relative importance of each component.

The techniques we developed are expected to perform well with sparse data sets as a consequence of previous findings on X-CNNs and the representations of time used, which permit potentially large and irregular gaps between data points. Our results are encouraging from both a computational perspective and a medical one: incorporating even simple descriptions of time into models can yield notable improvements, with more suitable temporal descriptions producing significantly better classification performance.

2 Medical Background

There has been growing interest in automated disease diagnosis and assessment based on machine learning (ML) models. This interest has focused on extracting features from time-series, such as using EEG data to predict epileptic attacks up to an hour in advance [6], and on various types of medical imaging. In [1], Ashinsky et al seek to predict osteoarthritis using images from MRI scans, whilst in [4], Cupek et al worked on the automated assessment of joint synovitis using ultrasound images, applying a set of image processing and ML techniques to determine synovitis severity. However, these two areas of interest have typically been studied separately, thereby potentially excluding relevant data, features, and modelling techniques which could improve model performance.

The case study explored in this paper concerns bone remodelling and the classification of imbalanced remodelling. Remodelling encompasses bone changes due to tissue maintenance and is governed by two interacting mechanisms: formation and resorption, in turn dependent on the RANK-RANKL-OPG signalling pathway [12]. Everyday actions, such as walking, cause micro-cracks in bone to develop. If this micro-damage were simply filled in with mineralised extracellular matrix, the old damaged tissue would continue to be weak, compromising the mechanical integrity of the bone and risking further micro-damage accumulation. Resorption is the process of removing the damaged tissue, triggered by osteocytes signalling and performed by osteoclasts. Osteoblasts prevent excessive bone excavation by producing a decoy receptor, OPG, as they mature, inhibiting osteoclastic activity. The osteoblasts then deposit collagen and induce mineralisation, thereafter differentiating themselves as osteocytes embedded in the extracellular matrix. In healthy bone tissue, these processes are coupled such that the site of the micro-crack is repaired completely. However, pathologies can cause imbalances that affect bone remodelling. Osteoarthritis and osteopetrosis cause excessive mineralisation, typically affecting joints due to these experiencing the greatest impacts and general wear-and-tear. This can ultimately lead to inflammation of synovial tissues, restricting movement and paining the sufferer. In contrast, diseases such as osteoporosis and osteomyelitis are caused by excessive osteoclastic activity gradually degrading bone density and reducing structural integrity, in turn making the damaged bone more vulnerable to further micro-cracks, causing additional deterioration in a vicious cycle. The standard approach to modelling the coupled processes of formation and resorption is to use stochastic simulations based on ordinary different equations or, when extending temporal models to spatio-temporal ones, partial differential equations. These are important for understanding systems and can use formal and probabilistic checking to verify simulation properties, but are limited by the complexity of extending such models and their inability to be used in automated data-processing tasks.

3 Mouse Data

The data were µCT scans of the tibiae of 15 mice taken over an eight-week period, between the pages of 14 and 22 weeks, which were used to study longitudinal bone changes caused by the processes of bone resorption and formation [8]. The 15 female mice were divided into three groups ($N=5$ per group), as shown in Figure 1, including a healthy control group and two groups with induced imbalanced bone remodelling. The control group were referred to by the researchers as being wild type. The bone formation group was treated with parathyroid hormone (PTH) — an anabolic peptide — four weeks into the experiment, at 18 weeks of age. As such, from ages 14 to 18 weeks this group was treated identically to the control group. The bone resorption group underwent ovariectomy at 14 weeks of age, causing œstrogen deficiency and accelerating bone resorption in the weeks 14 to 22, particularly localised in the most proximal part of the bone, close to the knee. Only female mice were considered for this reason, as observed differences could be attributed to treatments rather than potentially being the result of the sex of the mice in the different groups. For our case study, we considered data from only the wild type and PTH groups.

In order to monitor bone changes over space and time, weekly µCT scans of the right tibia of each mouse were taken. These high-resolution scans were reconstructed in a series of slices: cross-sectional images perpendicular to the axis of the tibia, as depicted in Figure 2, with a voxel size of 10µm in between the growth plates of the tibiae. This resulted in around 1200-1400 slices per mouse per week, giving a total of around 52 000 images with a combined size of around 20GiB per test group. Two images, both from the same mouse at the same time, but at different longitudinal positions, are given in Figure 3. Dr Dall’Ara’s research team performed image registration to place the images from different mice and different time points in the same reference system [11]. Due to the use of an operator crop to remove regions which did not include the tibia from the images generated by the µCT scanner, and in part due also to the mice still growing slightly, despite being considered skeletally mature, the image dimensions differ slightly from week to week and from mouse to mouse. This variation in dimensions is characterised by image heights between 400-500 pixels and widths between 500-700 pixels. All images were provided and stored in DICOM (Digital Imaging and Communications in Medicine) format. Many image processing libraries cannot accommodate variable image size, and the data were incomplete: no data were recorded for any of the mice for week two of the experiment, and one of the mice in the induced bone formation (anabolic) group was missing data for the final week. Furthermore, the data exhibited a issue common in medicine and bioinformatics: being deep rather than wide.

4 Preprocessing

In order to both accommodate the differing image sizes and week-to-week alignment variations, and to produce the different temporal descriptions, a number of pre-processing steps were used. The pipeline of these steps is shown in Figure 4.

As most image processing libraries require images with consistent dimensions, it was necessary to convert all images to have the same width and height. As the full dataset was known in advance, it was possible to standardise the dimensions by expanding each image to the maximum height and width values observed in the dataset, which were 501 and 763 pixels respectively. This expansion process placed the image centrally and set any surrounding pixels to the minimal intensity value in that image. When the maximal image dimensions are not known in advance, images could be up- or down-sampled to fit within pre-set dimensions, either preserving relative dimensions and filling empty pixels as described, or by fitting the image to the pre-set dimensions without preserving relative sizes. This resizing process was performed once and resulting images saved to disk to avoid the cost of resizing during each training phase. In order to perform differencing between any two expanded images, it was first necessary to align them, as the mouse from which the images were taken would not be in exactly the same position from week to week. The alignment was performed using scikit-image’s registration method, which for a 2D image outputs a pair of numbers indicating the $x$ and $y$ axis translations to perform. The first image underwent the specified translation for optimal alignment with the second image, followed by pixel-wise subtraction to create the difference map.

When each image was processed, an integer timestamp was recorded in the image meta-data, representing the week the image was taken in since the experiment start. Integer week numbers were recorded rather than normalised values for generalisability beyond the dataset of the case study. Timestamps normalised to the range [0,1], for example, this would require experiment durations to be known in advance, thereby excluding later data. Such timestamps would furthermore be incompatible between datasets which, whilst otherwise comparable, had been collected over different durations.

The final stage of preprocessing was to split the dataset into three partitions for training, validation, and testing. The data for one mouse of the five in each treatment group was selected at random and held out solely for testing. The remaining data was randomly permuted and partitioned as $90\%$ for training and $10\%$ for validation. Random permutation helps to avoid biases from repeated exposure of a network to similar training sequences and their associated gradient updates, thereby mitigating the risk of overfitting.

Due to the nature of the images, some common preprocessing steps were either undesirable or unnecessary. As the µCT scans were monochromatic, there was no opportunity to split the images into multiple colour channels or use intensity-chromacity encodings. Separating the differing information contained in these image channels can boost model performance [16], thus, where applicable, chromacity and/or intensity channels should be incorporated into medical X-CNNs, applied to both the original images and the difference images. Doing so would require only minor modifications to the X-CNN approaches presented herein. In order to prevent compression artifacts, no compression was applied to the images. Instead, the neural networks effectively perform sub-sampling via pooling, in the process learning weights; sub-sampling external to the network removes learning opportunities for the network itself. Finally, data augmentation - the process of synthesising novel datasets via some perturbation or transformation of an existing dataset - was not applicable. Augmentation is fundamentally incompatible with differencing approaches as the synthesised images are generated independently from one another. This means there is no guaranteed correlation in the differences between a reference-comparison image pair, so performing differencing on augmented images would amount to differencing randomly-selected image pairs, thereby invalidating the intent of differencing. Furthermore, augmentation is frequently used to expand small datasets and “fill in” sparse ones. For the dataset in this case study, there were large quantities of data for each mouse - around $10000$ images - and as the µCT process produced image slices with high spatial proximity, adjacent images should be highly correlated and with minimal variation. Effectively, adjacent image slices appear to be perturbed, noisy versions of one another and week-to-week alignment variations effectively introduce small translations and rotations, both whilst retaining the temporal relationships between image slices of the same individuals from different timestamps. Thus, augmentation would be unlikely to benefit this and similar datasets.

5 Models

We present three different temporal descriptions and two combinations thereof, compared against a state-of-the-art CNN baseline. Whilst the architectures for each approach differ, they share common elements. All models used a batch size of 32, employed max-pooling, batch normalisation, and dropout, and utilised rectified linear unit (ReLU) activation functions for their strong performance with convolutional networks [9]. We chose He initialisations for suitability with ReLU activations [7], as the Glorot initialisation scheme on which they are based, whilst generally applicable, was designed for sigmoidal activations [5]. The output layer of each network was identical: one neuron per treatment class, values normalised by softmax, with categorical cross-entropy as the loss function.

The basic building block of the convolutional part of each network is the chain, shown for non-cross modal networks in Figure 5 and for cross-modal networks in Figure 6. Chains consist of one or more convolutional layers, followed by pooling and any optional regularisation layers - in this case, batch normalisation and dropout. In the fully-connected layers, it was found that batch normalisation and dropout conflicted, with batch normalisation on its own providing the best results. Increasing the number of convolutional chains also proved to be unsuccessful, with classification performance dropping as more chains were added. This suggests the features of use to the network are subtle and that further compression of them into more compact features by pooling loses information. This supports the decision to not downsample the images during preprocessing. Furthermore, it was found that increasing the number of neurons in the first perceptron layer to 128 or 256 did not improve the models. This was in fact advantageous, as the majority of trainable parameters in each model,  $96\%$ for the differencing approaches, came from the connections between the convolutional and perceptron layers. Using fewer neurons in the first perceptron layer made the networks far smaller, allowing faster training.

In the case of X-CNNs, cross-connections are appended to the chain structure from Figure 5, creating what we refer to as X-chains, which allow the exchange of feature maps between channels. These cross-connections consist of $1\times 1$ convolutions of each channel, merged via tensor concatenation. Regularisation layers such as batch normalisation may happen either before or after the cross-connection process; the impact of this is not yet fully understood.

5.1 Baseline CNN

The baseline we used for comparison was a CNN, incorporating state-of-the-art techniques. There were five convolutional chains, the last of which did not use dropout and fed to a series of fully-connected, “dense” layers, comprising 64, 32, then $n$ neurons, where $n$ was the number of treatment classes – 2 in this case study. All convolutions used $3\times 3$ kernels – 16 such kernels for the first chain and 32 for each subsequent one. The first two fully-connected layers were regularised by batch normalisation on their outputs. The max-pooling layers in the chains all used $2\times 2$ non-overlapping blocks, the dropout probability was $0.5$, and $\gamma=0.003$ for L2 regularisation in the fully-connected layers. As mentioned, kernel initialisations for all but the output layer were He Gaussians, accompanied by ReLU activations.

5.2 Timestamps

The most obvious and simple approach to incorporating time is the use of a single number indicating the timestamp, which may be cardinal or ordinal. As has been discussed, timestamps were not normalised. The timestamp is provided via a unit-sized tensor, concatenated with the output of the convolutional layers as input to the first fully-connected layer. Otherwise, this model was identical to the baseline.

5.3 Difference Images

Difference images provide an implicit description of time by performing a pixel-wise subtraction between an image at the time step, $t$, of interest — the comparison image — and an image at a preceding time step $t_{0}$ — the reference image, should one exist:

[TABLE]

Two approaches to selecting the reference image are to use images from the first recorded timestamp for each mouse or individual, and to use images from the most recent timestamp available for each mouse or individual. These are referred to as absolute and relative differencing respectively, as the timestamps selected are either absolute or are relative to the current timestamp. For either approach, the first timestamp’s difference images will always be of uniform intensity equal to a minimum value; this represents the lack of preceding data. That minimum value may be known in advance, e.g. [math], or may be the minimum value observed in each comparison image. The former requires pre-calculation or determination of a minimum value, which may not always be possible, and may introduce values in the difference image much larger than would be obtained in later timestamps’ difference images, placing undue emphasis on the first timestamp and the least meaningful difference image. Furthermore, it is inconsistent with the padding approach used in the image expansion process in section 4. The latter approach, however, introduces global variation by allowing different intitial-timestamp difference images to have different values. We selected the latter approach as it is more suitable for large and real-time data sets. Differencing faced one further difficulty, in that for each mouse there were a different number of images each week. As such, a one-to-one matching between reference and comparison images was not possible. As each µCT scan began at the proximal growth plate, we paired the $k^{th}$ reference image with the $k^{th}$ comparison image from the appropriate week, stopping when there were no more images from either of the weeks. This lead to $2.18\%$ of images being unused.

Having both original and difference images provides multiple processing channels thereby enabling the use of cross-connections, making such models X-CNNs. Processing was identical for absolute and relative differencing so that any observed difference in performance could be attributed to the type of differencing used, rather than potentially being the result of different architectures or hyperparameters. The processing in both convolutional channels used the same structure and hyperparameters, namely five X-chains as in Figure 6 with dropout probability set to $0.25$, He Gaussian kernel initialisations with ReLU activations, and max-pooling over non-overlapping blocks with stride length $2$. The number of convolution kernels used in each channel was eight for the first two chains, 16 for the subsequent two chains, and 32 for the final one. The first two fully-connected perceptron layers contained 64 and 32 neurons respectively, both using L2 regularisation with $\lambda=0.0003$. The final, output, perceptron layer was identical to that used in the other models.

5.4 Combining Descriptions of Time

It is, of course, possible to combine both explicit descriptions of time, like timestamps, with implicit ones, like differencing. Providing a network with both representations allows it to determine the relative importance of each during training, thereby adapting better to the particular domain and problem. There may be situations where one temporal descriptor is more useful than another or when one can be a deciding factor.

For combined temporal descriptions, we have three separate input channels: the original image, the difference image, and the timestamp. The first two undergo the same processing as in section 5.3, then the tensors from all three channels are concatenated before input to the fully-connected perceptron layers. This is shown in Figure 7. All hyperparameter choices were as in section 5.3.

6 Results

In Table 1, we present the test-set accuracies and mean training times over five and 10 epochs, together with the number of parameters (both trainable and non-trainable), for each of the six models considered. We use accuracy as an evaluation metric as the training, validation, and test sets were all almost perfectly balanced in terms of the number of images available for each treatment class. The exact splits are shown in Table 2. We provide the average training time per epoch, in addition to the number of epochs of training, as real-world training times for new models may impact decisions on which models to use in a medical setting based on available resources and how frequently models need to be re-trained.

The best performing models were those using absolute differencing, achieving over $95\%$ accuracy after both 5 and 10 epochs of training, both with and without timestamps. The inclusion of timestamps with absolute differencing yielded the best results overall - $99.54\%$ accuracy after 10 epochs. In contrast, the CNN baseline managed only $73.02\%$ accuracy, achieved with 5 epochs of training. Whilst not out of line with the performance of many classification and prediction tasks on medical datasets [1, 6, 13, 14], the relatively poor performance of the baseline when compared to commonplace results on datasets such as MNIST and CIFAR-10 highlights the inherent difficulty of medical tasks and the necessity of finding novel, informative features. We see that, for the same number of training epochs, spatio-temporal models always outperformed the purely spatial baseline and that only in one of the 10 results for the temporal models did the accuracy fall below a baseline result: the CNN with timestamps, trained for 10 epochs, performed worse than the baseline CNN trained for five epochs. In general, it was the case that training for longer, 10 rather than five epochs, offered little to no benefit — for only two of the six models was there a performance increase, and then only of circa $1\%$, whereas the other models saw drops of around $3-11\%$, averaging a $-7.37\%$ decrease. This suggests the models may be susceptible to overfitting, at least on a deep rather than wide dataset as in the case study, despite extensive use of regularisation techniques. Given that performance after only five epochs could reach around $98\%$, this is in fact positive – powerful, useful models can be trained more quickly. Note that if difference images provided no information beyond that contained in the original images, we would expect to find that training the non-differencing models for twice as many epochs would allow those models to reach the same accuracies as with differencing, as they would have processed the same total quantity of information. As this was not the case, we must conclude that difference images contain valuable additional information.

Whilst relative differencing exceeded the baseline and timestamp accuracies, it performed noticeably worse than absolute differencing — between approximately $10-25\%$ worse, depending on the number of training epochs and the inclusion or exclusion of timestamps. This is likely a consequence of the dependency range the model can represent, with relative differencing only capturing short-range dependencies in a limited temporal region, wheras absolute differencing is capable of modelling arbitrarily long-range dependencies as RNNs can, although more efficiently as RNNs build long-range dependencies from iterating over short-range ones. For conditions like abnormal bone remodelling which, when averaged over cycles of periodic behaviour, exhibit monotonic trends in a metric of interest, such as bone density or volume, long-range dependencies may best capture these trends. Relative differencing may be more suitable for conditions characterised by fluctuations around a central (range of) value(s), which may average out over longer time spans. The efficacy of timestamps may be related to this, in that for absolute differencing, timestamps reinforce the information provided by the difference images - larger pixel values, indicating larger differences in bone density, tend to correlate with larger timestamps. On the other hand, relative differencing may produce similar difference images for each week, thereby making it difficult to distinguish which week a difference image was from; timestamps would not correlate with the pixel intensities in the difference images, thus attempting to coordinate the information from these two temporal descriptions may produce erroneous associations. This would explain the impaired performance when using relative differencing with timestamps.

From Table 1, we see that the non-cross modal models required fewer than 900K parameters each, whilst the cross-modal ones required a little over 1.5 million each. In comparison with recent high-performance image classification networks — 5 million parameters in GoogLeNet, 23 million in Inception v3, 60 million in AlexNet, and 180 million in VGGNet [15] — a 1–1.5 million parameter network is computationally inexpensive. Based on the training times in Table 1, we see two main points. Firstly, that the difference image-based approaches take around four times as long as the non-differencing ones, likely due to a combination of disk-load speeds from a larger-than-memory dataset and the cost of training a more complex network with around twice as many parameters. Secondly, training these networks necessitates the use of GPUs or highly parallel multi-core processing systems — whilst smaller than many high-performance models, they are not sufficiently computationally inexpensive on train on standard desktop systems.

7 Conclusions

We have shown several approaches to the incorporation of temporal information for image-processing tasks and compared their performance on a balanced dataset comprising over 100 000 data points across two treatment groups. For the same number of epochs of training, the models with temporal descriptions always outperformed an atemporal state-of-the-art CNN baseline, in some cases by in excess of $25\%$, highlighting the distinct potential to improve atemporal models of temporally-varying processes and systems. Whilst even simple mechanisms for providing temporal information can yield benefits, such as integer timestamps, the best performance is achieved by using a temporal description suitable for the domain and problem. For our case study, wherein the treatments produced monotonic effects over extended periods of time, the ability of absolute differencing to capture long-range dependencies produced the best results. For certain medical applications this is highly advantageous — accurately determining the disease or treatment status of a patient may require as little as two scans some number of weeks apart, without the need for invasive procedures or regular appointments. Furthermore, the temporal models we present permit sparse, non-sequential data sampled at irregular intervals, making them far more flexible than RNNs and 3D CNNs and without the need for interpolation or other approaches to reconstruct sequential data, which ultimately add a preprocessing cost for those types of networks.

Whilst we have shown the effectiveness of cross-modal neural networks for a particular temporally-varying medical dataset, our approach was kept general throughout and is applicable to a wider medical context. The models we developed could easily be extended or modified to particular domains or for further research into alternative temporal descriptions. Potential extensions include the use of both absolute and relative differencing, to capture both long- and short-range dependencies, the use of SVMs or RNNs either in place of or subsequent to the fully-connected layers, and the application of our models to colour images. The latter would allow investigation of whether differencing extracts additional information from colour channels or only from the original image or luminance channel. In order to create less computationally-expensive and resource-intensive networks, the use of only difference images without their corresponding reference images should be explored, as they retain some amount of spatial information whilst representing temporal variations.

The proposed approach has the potential for improving the current assessment of the effects of treatments for musculoskeletal pathologies preclinically. The classification of the significant effects of a treatment can be used to optimise the dosage, time of treatment, or even combined interventions, for example alternate or overlapped treatment with anabolic and anti-resorptive drugs, or usage in combination with mechanical stimulation. Longitudinal experiments in mice are expensive and have ethical issues related to the usage of animals in research. Therefore, improving the classification of the effects of interventions with an automated, operator-insensitive tool would lead to a dramatic reduction in cost and the number of animals in research, in line with the 3 Rs: reduction, refinement, and replacement.

From the perspective of the longitudinal characterisation of phenotypes, we believe that there are two main challenges for the future. Firstly, if researchers wish to have broad oversight on the effects of ageing, co-mordibities, and their related interventions rather than focusing on specific tissues, organs, or systems, then it will be vital to combine data obtained from different modalities. As the data collected from different modalities is unlikely to be independent, uncovering powerful, effective ways to combine this information should result in models which are greater than the sum of their parts. Differencing goes some way towards this by using a spatial format — an image — to represent a change over time, and cross-modal networks support feature-sharing between selected modalities, allowing researchers to encode domain-specific knowledge via carefully chosen cross-connections. The second challenge we perceive is the need for increased image resolution to capture features which are not presently available, such as osteocyte lacunae within the extracellular matrix. The availability of such high-resolution images and the features they would unlock would require research on how to create machine learning approaches capable of accommodating the width, variety, and detail of such features in an efficient, flexible manner, with consideration for the likelihood of data points being sampled irregularly, at different rates, and with missing values.