A Convolutional Neural Network for the Automatic Diagnosis of Collagen VI related Muscular Dystrophies

TL;DR

This paper presents a CNN-based computer-aided diagnosis system for rare Collagen VI muscular dystrophies using confocal microscopy images, addressing data scarcity and aiding diagnosis and therapy monitoring.

Contribution

It introduces a patch-based CNN approach for diagnosing low-prevalence muscular dystrophies from microscopy images, providing localized and global diagnostic insights.

Findings

Effective patch classification for disease detection

Provides both localized and overall diagnostic decisions

Supports monitoring therapy progress

Abstract

The development of machine learning systems for the diagnosis of rare diseases is challenging mainly due the lack of data to study them. Despite this challenge, this paper proposes a system for the Computer Aided Diagnosis (CAD) of low-prevalence, congenital muscular dystrophies from confocal microscopy images. The proposed CAD system relies on a Convolutional Neural Network (CNN) which performs an independent classification for non-overlapping patches tiling the input image, and generates an overall decision summarizing the individual decisions for the patches on the query image. This decision scheme points to the possibly problematic areas in the input images and provides a global quantitative evaluation of the state of the patients, which is fundamental for diagnosis and to monitor the efficiency of therapies.

Tables

Table 1. Table 1: The details of the CNN architecture proposed in this paper.

Layer	Type	Number of neurons (output layer)	Kernel size	Stride
1	Convolution	64x64x128	3x3	1
2	Max Pooling	32x32x128	2x2	2
3	Convolution	32x32x64	3x3	1
4	Max Pooling	16x16x64	2x2	2
5	Convolution	16x16x32	3x3	1
6	Max Pooling	8x8x32	2x2	2
7	Fully Connected	150	-	-
8	Fully Connected	2	-	-

Table 2. Table 2: The performance of the proposed systems on the test set.

System	Size	$A$	$P$	$R$	$F_1$
Patch	64x64	0.95	0.92	0.98	0.95
Image	1024x1024	1	1	1	1

Equations

$$y=\sum _{n=1}^N{I}_n{k}_n,$$

$$f(x)=\left\{\begin{array}[]{ll}0&x<0\\ x&x\geq 0\end{array}\right..$$

$$\sigma (x)=\frac{1}{1+e^{-x}},$$

$$E=\frac{1}{2}\sum _{i=1}^n||o_i(W)-t_i|{|}^2,$$

$$W=W-\alpha \frac{\partial E}{\partial W},$$

$$A=\frac{t_p+t_n}{t_p+f_p+t_n+f_n},$$

$$P=\frac{t_p}{t_p+f_p},$$

$$R=\frac{t_p}{t_p+f_n},$$

$$F_1=\frac{2PR}{P+R},$$

Full text

††journal: Expert Systems With Applications

1 Introduction

Deficiencies in the structure of collagen VI are a common cause of neuromuscular diseases with manifestations ranging from the Bethlem myopathy to the severe Ullrich congenital muscular dystrophy. The symptoms include proximal and axial muscle weakness, distal hyperlaxity, joint contractures, and critical respiratory insufficiency, which requires assisted ventilation and results in a reduced live expectancy. Moreover, the skin and other connective tissues where collagen VI is abundant are also affected [1, 2]. The collagen VI structural defects are related to mutations of three main genes (COL6A1, COL6A2, and COL6A3, OMIM 254090 and 158810). Thus, the new advances in genome editing tools open the possibility to successfully treat these neuromuscular diseases for the first time. This opportunity, though, comes with important challenges. Beyond the challenges of gene editing, this paper focuses on the challenges arising when trying to formally evaluate the efficiency of the therapeutic approaches in the recovery of the collagen VI microfibrillar network.

Despite the implementation of next generation of gene sequencing, the diagnosis of rare diseases remains challenging. This is particularly true for diseases caused by dominant mutations where there is not a complete absence of a main protein and when the effect of a genetic variant on the protein structure may not be evident. Thus, prior to genetic analysis or to demonstrate the pathogenicity of novel mutations, the standard technique for the diagnosis of collagen VI related dystrophies is the analysis of images of fibroblast cultures [3] (see Fig. 1). Several aspects of the images such as the orientation of the collagen fibers, the regularity of the collagen network, and the arrangement of cells in such network are taken into account by the specialists to identify the patients. However, this evaluation is only qualitative and the regulatory agencies would not approve any treatment without an objective proof of its effectiveness [4]. Thus, there is an acute need of accurate methodologies to quantitatively monitor the effects of any possible new therapy.

This paper proposes a computer aided diagnosis (CAD) system for dystrophies caused by defects in the structure of collagen VI. The system uses a Convolutional Neural Network (CNN) to classify fibroblast culture images from possible patients. The CNN is trained from cases labeled by specialists but, as in the case of all rare diseases, there is a severe lack of training data. To address this issue, this paper uses a patch-based data-augmentation technique similar to the one proposed in [5].

This paper is structured as follows: Section 2 frames the contribution of this paper in the context of CAD systems using medical images and Section 3 introduces the basic elements defining a CNN. Then, Section 4 describes in detail the proposed system, Section 5 discusses the image acquisition and the data augmentation scheme, and Section 6 presents the results. Finally, Section 7 concludes the paper and points to issues deserving further attention.

2 Related work

The automatic analysis of medical images is typically divided in two steps. In the first one, relevant features are extracted from the image. In the second step, those features are used to evaluate a particular function that in our case is the classification between control and patient images.

Several methodologies have been proposed to identify relevant features in medical images. For instance, in [6] a multi-scale Gaussian filterbank is used to find perturbations in chest X-ray images. However, in our case, the disease causes a global disorganization of the collagen network, rather than local perturbations on it. In [7] the authors use a principal component analysis technique to extract relevant features from an X-ray image. These features, however are linearly related to the input image and such a simple transformation is not likely to provide relevant features in our case. In other approaches, Localy Binary Patterns (LBP) [8], completed LBP [9], gray-level co-occurence matrices [10], or Artificial Neural Networks (ANNs) [11] have been also used to obtain features, all on input images significantly different from the ones available in our case. In works more related to the problem considered in this paper, some systems have been proposed to extract features from collagen images, such as the distribution of the fibers orientation [12]. However, in our case, this distribution is only weakly related with the diagnostic, since some patients present coherently-oriented collagen networks (see Fig. 1-center). The same weak relation occurs with other possible features such as the average intensity of the collagen network, or the presence of collagen precipitates in the images. Thus, in the addressed problem, the selection of relevant features for the diagnostic is a challenge on itself.

Nevertheless, assuming that relevant features can be somehow identified, a classical approach would apply a classification method relying on them. Since there is not a single method adequate for all problems, different classifiers have been used on features extracted from medical images, including naive Bayes classifiers for the detection of brain tumors [13], support vector machines for the classification of pigmented skin lesions [14] and breast cancer [15], random forests for the diagnosis of tuberculosis [16], non-parametric approaches to identify patients with diabetic retinopathy pathology [17], and Convolutional Neural Networks (CNNs) for the detection of melanoma or glaucoma [18, 19]. In [20] some of these approaches are compared on histopathological breast cancer images, which are similar to the images considered in this paper, and they conclude that non-parametric approaches and CNNs provide the best results. CNNs are an evolution of the classical ANNs. However, there exists two main differences between ANNs and CNNs. First, CNNs are invariant to translations of the features in the image due to shared parameters between the image convolutions. Second, new types of layers appear in the CNNs such as the so called max-pooling layers, which enlarge the receptive fields of the convolutions and also contribute to the translation invariance. Such improvements, together with the availability of specialized hardware and software for training CNNs have established them as the de-facto standard tool for image processing, provided that enough training data is available. An additional advantage of the CNNs is that they can directly operate on images, taking care of both identifying relevant features and performing the classification relying on them [21, 22]. This is the fundamental reason to apply them to the problem addressed in this paper.

3 Preliminaries on CNNs

A CNN is composed by a set of layers of simple computational elements, called neurons, connected to elements in the previous layers. The layers are evaluated in sequence, from the input to the output. The computation elements are pre-defined and fixed and, thus, the parameters encoding the learned function are the weights associated with the connections in the network.

Three different kinds of layers are typically present in a CNN: convolutional layers, pooling layers, and fully connected layers [23]. Convolutional layers extract features applying a convolutional kernel all over the input image. Formally, each kernel is computed as

[TABLE]

where $I_{n}$ is a pixel in the image, $k_{n}$ is $n$-th kernel weight, and $N$ is the kernel size, i.e., the amount of pixels in the image affected by the kernel. Convolutions are typically applied on squared areas around selected pixels in the input image and, thus, the result of this operation, is a new image, $Y$, which is called a feature map. Often, the output of a kernel is passed through an activation function. For instance, in our approach, each kernel is followed by a rectified linear unit (ReLU) activation function,

[TABLE]

which eliminates negative inputs, introducing non-linearities in the CNN with a low computational cost. The horizontal/vertical distance (in pixels) between the centers of two consecutive kernels is known as the stride of the kernel. Strides larger than one are typically used to downsample a feature map. Alternatively, the pooling layers down-sample the input by summarizing a patch in the image with a single value. For instance, in our work we use the max-pooling operation, where the output is the maximum of the inputs.

Finally, the fully-connected layers are typically placed at the end of the network and they provide the classification decision by applying an activation function. In our case, we use the sigmoid activation function

[TABLE]

where $x$ is a weighted sum of inputs and where the output is a number between 0 and 1.

As depicted in Fig. 2, a CNN learns comparing its outputs with the ones given in a labeled dataset. Specifically, we want to make the output $o_{i}$ and the target $t_{i}$ to be as similar as possible for each input instance $i$ = 1,…,$n$, where $n$ is the number of instances. More formally, we want to minimize the error

[TABLE]

with $W$ the weights in the CNN determining the output. The weights are initialized with random values and they are updated with a gradient descent rule as

[TABLE]

where $\alpha$ is the learning rate determining the intensity of the gradient descent during training. The backpropagation algorithm [24] is an efficient way to compute the gradient of the error with respect to the weights and it is the approach used in this work.

For binary classification tasks, the available training data is split in two sets: a training set typically containing about 80% of the data and a test set. A particular batch of training inputs is used at each iteration and the process is repeated until all the training inputs are used. The iteration over all the available training set is known as an epoch. After each epoch, the performance of the CNN is evaluated using a validation set including 10% of the training data. Finally, after training for a certain number of epochs, the CNN generalization capability is evaluated using the testing set, that was not used for training.

The evaluation with a validation and a testing set reduces overfitting, that is a major issue happening when a CNN perfectly classifies the training data, but is unable to correctly classify unseen cases. Dropout is another mechanism to prevent overfitting [25]. With this mechanism, some neurons are deactivated with a given probability. These elements do not participate in the activation of posterior neurons nor in the error correction, i.e., they are virtually removed from the CNN. The benefit is that the network becomes less sensitive to the specific weights.

4 The proposed system

Figure 3 provides a global view of the proposed CAD system. The system is divided in four modules. The first module receives a full image and splits it into non-overlapping patches of 64x64 pixels. The second module is formed by the CNN classification model, that receives the patches and outputs an independent prediction for each one of them. The third module receives the local decisions for each patch and takes a global decision using majority voting. The last module visualizes the input image, the decision for each patch represented in a color code, and the overall decision of the system (see Fig. 4). Cyan is used to frame patches with more than 90% probability of belonging to the control class, steel blue is used for patches with probability between 70% and 90%, yellow for patches with probability between 50% and 70%, orange for patches with probability between 30% and 50%, and finally red for patches with less than 30% of probability of belonging to the control class. This color code offers the possibility of easily spotting suspicious areas in the image. The system also provides the overall decision on the image and a global score computed as the percentage of patches classified as control in the image. This score enables to track the evolution of a patient and, thus, it offers a tool to assess whether or not a particular treatment is effective.

The core component in the proposed approach is the CNN model used to classify the image patches. We experimentally identified a CNN architecture suitable for the addressed problem. In this experimental selection two concepts were considered: the size of the network (with small/large number of convolution kernels) and the abstraction level (with increasing/decreasing sizes of the subsequent layers of the CNN). Experimentally, small CNNs with decreasing number of features provided the best results. A key feature in such architectures, is that many low-level features are learned on the first layers and few high-level discriminating features of the images are generated in deeper layers. The proposed architecture is detailed in Table 1 and is similar to the one in [5], which has already been proven to be particularly adequate for image classification. However, our network is smaller since the classification task addressed here is simpler than the one addressed in [5].

Figure 5 provides a detailed overview of the CNN architecture proposed in this paper. The input image patch is normalized in order to have zero mean and unit variance. Normalization improves the gradient descent process and avoids its premature convergence. The normalized image patch is passed through three convolutional layers (layers 1, 3, and 5) which include the application of ReLU activation functions after the kernel computation. The first convolution layer defines 128 feature maps of size 64x64, the second one defines 64 feature maps of size 32x32, and the last one defines 32 feature maps of size 16x16. The reduction in the size of the feature maps is obtained with max-pooling layers (layers 2, 4, and 6) with a stride of 2 following the convolution layers. After the feature generation layers, the classification is implemented with two fully-connected layers (layers 7 and 8). The first one has 150 neurons, also followed by a ReLU activation function trained with a dropout mechanism with probability 0.5. The second fully connected layer has 2 neurons whose output is truncated into a single binary output by a sigmoid activation function to provide the final classification.

5 Data acquisition and augmentation

Samples from the forearm were obtained from possible patients with as well as from aged-matched controls. Primary fibroblasts cultures were established using standard procedures [3]. Confluent fibroblasts (patient and control in parallel) were treated with 25 $\mu$g/mL of L-ascorbic acid phosphate magnesium (Wako Chemicals GmbH, Neuss, Germany) for 24 hours. After that time, cells were fixed with 4% paraformaldehyde in phosphate-buffered saline solution. Collagen VI was detected by indirect immunofluorescence using a monoclonal antibody (MAB1944, Merck, Germany) as previously described in [3] and fibroblast nuclei were stained using 4,6-diamidino-2-phenylindole (Sigma Chemical, St. Louis, USA).

The images were acquired with a Leica TCS SP8 X White Light Laser confocal microscope with hybrid spectral detectors (Leica Microsystems, Wetzlar, Germany). The confocal images were acquired using a HC x PL APO 20x/0.75 dry objective and with the confocal pinhole set to 1 Airy unit. Collagen VI was excited with an argon laser (488 nm) and detected in the 500-560 nm and nuclei were excited with a blue diode laser (405 nm) and detected in the 420–460 nm. Appropriate negative controls were used to adjust confocal settings to avoid non-specific fluorescence artifacts. The detector gain and offset values were adjusted to use the entire dynamic rate of the detector (12 bits) and to avoid oversaturated voxels. Sequential acquisition settings were used to avoid inter-channel cross-talk. Ten sections of each sample were acquired every 1.5 $\mu$m along the focal axis (Z-stack) and combined into an integrated intensity projection to form a single image.

Since the data acquisition procedure is complex and we are dealing with a rare disease, the available data is limited. Thus, we propose a data augmentation scheme inspired in the one described in [5] to generate enough inputs to train the CNN.

The data is augmented by splitting the given full size images into small, non-overlapping patches. Each patch is used as an independent input to train the CNN model. In our approach, the patch is of size 64x64 pixels since it has been shown to be a particularly relevant window size for CNN classification models on similar images [26]. In our case, the patches capture the main features of the input images, i.e., a significant portion of the collagen network and several nuclei or parts of them.

Since our input images are of 1024x1024 pixels, 256 patches per image are generated. Each patch is further transformed to get even more variations of the data: They are rotated clockwise by 90, 180, 270 and 360 degrees and every rotated patch is flipped horizontally. Thus, eight different variations of each original patch can be obtained and, consequently, each input image generates 2048 training inputs for the CNN. Since initially we have 276 images, with this data augmentation process, the training and testing sets include, respectively, 56320 and 14336 patches, without taking into account the rotated and mirrored patches.

6 Results and discussion

The proposed system was implemented using the Python programming language and the Keras 2.2 library with TensorFlow 0.19 backend and it was executed on a workstation with a Intel Core i7-7700HQ processor and 16 GB of RAM. All the processing was performed in a NVIDIA GeForce GTX 1050 GPU. The batch size is set to 32 and the training is executed while the performance of the classification of the test set improves. On average, each epoch took about 14 minutes to complete and less than 10 epochs are necessary to converge. For training the neural network we use the Adam [27] optimizer, with a learning rate of $\alpha=0.01$, an exponential decay rate for the first moment estimates of $\beta_{1}=0.9$, an exponential decay rate for the second moment estimates of $\beta_{2}=0.999$, and $\epsilon=10^{-8}$ to prevent any division by zero.

We present results at two different levels. The first level is the per-patch performance of the CNN model and the second is the performance of the proposed system on full images, by using majority voting. In this way, we evaluate how the model performs using only local features (patches) and how it improves when integrating local results to take a global decision (full images).

Figure 6 gives the confusion matrix of the system diagnosing on 64x64 patches on the test set. In this confusion matrix, the amount of patches correctly and incorrectly classified as belonging to the control class are in the first row. From a total of 8115 control patches in the test set, 7585 are true negatives ($t_{n}$), i.e., the inputs correctly classified as control and 530 are false positives ($f_{p}$), i.e., the inputs incorrectly classified as patient. The classification of patients is given in the second row, where, 102 are false negatives ($f_{n}$), i.e., the inputs incorrectly classified as control, and 6119 are true positives ($t_{p}$), i.e., the inputs correctly classified as patient. In the majority voting system the accuracy is perfect and, thus, the confusion matrix is trivial and not given here.

Table 2 gives the performance evaluation metrics for each of the analyzed levels. The accuracy, $A$, precision, $P$, recall, $R$, and $F_{1}$ score are used to assess the classification performance. The accuracy,

[TABLE]

refers to the correct classification rate, defined as the ratio of correctly classified cases with respect to the total number of cases. The precision,

[TABLE]

is the ratio of inputs correctly classified as patient with respect to the total number of inputs (correctly or incorrectly) classified a patient. The recall,

[TABLE]

gives the ratio of samples correctly classified as patient with respect to total inputs that are actually in the patient class. Finally, the $F_{1}$ metric

[TABLE]

is the weighted harmonic mean of precision and recall.

At the level of patch classification, the system achieves an accuracy of 0.95, a precision of 0.92, and a recall of 0.98. Consequently, the proposed system attains a $F_{1}$ score of 0.95. At the level of images, the system achieves perfect classification results.

CNNs can provide very accurate results, but they are often seen as black boxes since it is hard to find out what are the features and the classification criteria learned from the data. However, recent works try to interpret the CNNs after the training phase [28]. Following this approach, we provide a visualization of two relevant parts of the CNN. First, Fig. 7 shows the contribution of each pixel in the image to the classification as control class. The network focus on the fibroblast nuclei and on the regions without collagen VI. Note that the intensity on the later is lower, but their extension is larger and, thus its actual contribution to the classification is also significant. In second term, Fig. 8 shows the saliency maps of the convolutional layers with respect to the patient class, where the saliency is the gradient of the output with respect to the corresponding convolutional layer. At early layers, the network learns to detect low level features, such as the presence of fibroblasts or the lack of collagen VI, confirming the previous analysis. In posterior layers, the features become more abstract and, thus, more difficult to relate with elements in the input images.

7 Conclusions

This paper describes a system for the computer aided diagnosis of muscular dystrophies caused by deficiencies in the structure of collagen VI. The proposed system relies on a deep convolutional neural network and on a data augmentation scheme to handle the problem of lack of data typical of rare diseases. The proposed system is capable of achieving perfect results in the diagnosis task in the available testing set. This provides a solid tool to automatically and accurately track the effect of potential therapies for the recovery of collagen VI in patients. Furthermore, we visualized the contribution of each pixel to the classification as well as the features learned by the CNN in each convolutional layer by means of saliency maps. This identifies the primary features discovered by the CNN. Our current research endeavours focus on extending this analysis to clarify the abstract features and the classification criteria identified by the proposed CNN. Moreover, we are studying the possible application of image-based computer aided diagnosis procedures to other rare diseases related to deficiencies in collagen VI or other extra-cellular matrix proteins.

Ethics Statement

This study was carried out in accordance with the recommendations of the Fundación Sant Joan de Déu Ethics Committee. Written informed consent was obtained from patients and/or their parents or guardians in accordance with the Declaration of Helsinki. The protocol was approved by the Fundación Sant Joan de Déu Ethics Committee.

Acknowledgements

Adrián Bazaga was supported by a JAE-intro scholarship granted by the Spanish Council of Scientific Research. Josep M. Porta is funded by the Spanish Ministry of Economy and Competitiveness under project DPI2017-88282-P.

Cecilia Jiménez-Mallebrera and Carmen Badosa are funded by the Health Institute ’Carlos III’ (ISCIII, Spain) and the European Regional Development Fund (ERDF/FEDER), ’A way of making Europe’, grants references PI16/00579, CP09/00011 and Fundación Noelia.

Author Contributions

Adrián Bazaga – Conceptualization; Investigation; Formal analysis; Methodology; Software; Validation; Visualization; Writing - original draft; Writing - review & editing.

Mònica Roldán – Data curation; Resources; Writing - original draft; Writing - review & editing.

Carmen Badosa – Data curation; Writing - review & editing.

Cecilia Jiménez-Mallebrera – Data curation; Resources; Writing - original draft; Writing - review & editing.

Josep M. Porta – Conceptualization; Formal analysis; Methodology; Supervision; Funding acquisition; Writing - original draft; Writing - review & editing.