Science, Faculty of
Zoology, Department of
DSpace
UBCV
Holmes, Nathan
Master of Science - MSc
University of British Columbia
Assessing the morphology of myelin in transmission electron micrograph (TEM) images is an important research tool for studying certain neuropathologies and neuroregeneration. Additionally, TEM can be used to validate the accuracy of other methods of measuring the state of myelin as it allows direct visualization of myelin ultra-structure. Often this requires laborious manual quantification of TEM micrographs to extract useful morphological data, especially in cases where the feature of interest is not commonly measured such that software is not commercially available to aid in quantification. A recent study produced by our lab involved comparing measurements of myelin based on magnetic resonance imaging (MRI) T2 relaxation with TEM images in conditions of healthy and degenerating myelin in the rat spinal cord following a dorsal column transection injury. The MRI T2 techniques can quantify the amount of water trapped between layers of the myelin sheaths, as opposed to water in other spaces in healthy tissue. However, what the MRI T2 measurement was actually corresponding to in degenerating tissue was unknown. The TEM portion of this analysis involved manually tracing the area of the myelin sheath, and thresholding the image to separate intact myelin from vacuous spaces occurring in the myelin as it degenerates. This made it possible to determine that the MRI T2 technique only measures to the quantity of intact myelin in degenerating tissue. This thesis uses the TEM images and the manually quantified intact myelin measurements as data-sets for training machine learning approaches to automate this task. Three approaches are trained and evaluated: a k-nearest neighbour (KNN) approach based on natural clustering of image histograms, a standard SegNet deep-learning artificial neural network, and a more generic dense prediction artificial neural network featuring skip-connections to preserve low level data. The KNN approach was found to be completely unsuitable for the task. The SegNet approach had substantial difficulties which indicate that it did not have the flexibility to fully learn the data-set. Finally, the generic dense prediction network performed robustly in both the training and test data-sets. Additional testing is performed to contrast the performance of the two artificial neural networks.
Application of Machine Learning Techniques for the Quantification of Myelin Ultra-structure inTransmission Electron Micrographs of the Rat Spinal CordbyNathan HolmesB.S., The University of British Columbia, 2011B.C.S., The University of British Columbia, 2017A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMaster of Sciencein THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Zoology)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)November 2023© Nathan Holmes, 2023The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled:Application of Machine Learning Techniques for the Quantification of Myelin Ultra-structure inTransmission Electron Micrographs of the Rat Spinal Cordsubmitted by Nathan Holmes in partial fulfilment of the requirements forthe degree of Master of Sciencein ZoologyExamining Committee:Wolfram Tetzlaff, Professor, Zoology, UBCSupervisor Piotr Kozlowski, Professor, Radiology, UBCSupervisory Committee Member,Doug Altshuler, Professor, Zoology, UBCAdditional ExaminerAdditional Supervisory Committee Members:Vanessa Auld, Professor, Zoology, UBCSupervisory Committee MemberiiAbstractAssessing the morphology of myelin in transmission electron micrograph (TEM) images is an important research tool for studying certain neuropathologies and neuroregeneration. Additionally, TEM can be used to validate the accuracy of other methods of measuring the state of myelin as it allows direct visualization of myelin ultra-structure. Often this requires laborious manual quantificationof TEM micrographs to extract useful morphological data, especially in cases where the feature of interest is not commonly measured such that software is not commercially available to aid in quantification.A recent study produced by our lab involved comparing measurements of myelin based on magnetic resonance imaging (MRI) T2 relaxation with TEM images in conditions of healthy and degenerating myelin in the rat spinal cord following a dorsal column transection injury. The MRI T2 techniques can quantify the amount of water trapped between layers of the myelin sheaths, as opposed to water in other spaces in healthy tissue. However, what the MRI T2 measurement was actually corresponding to in degenerating tissue was unknown. The TEM portion of this analysis involved manually tracing the area of the myelin sheath, and thresholding the image to separate intact myelin from vacuous spaces occurring in the myelin as it degenerates. This made it possible to determine that the MRI T2 techniqueonly measures to the quantity of intact myelin in degenerating tissue. This thesis uses the TEM images and the manually quantified intact myelin measurements as data-sets for training machine learning approaches to automate this task. Three approaches are trained and evaluated: a k-nearest neighbour (KNN) approach based on natural clustering of image histograms, a standard SegNet deep-learning artificial neural network, and a more generic dense prediction artificial neural network featuring skip-connections to preserve low level data. The KNN approach was found to be completely unsuitable for the task. The SegNet approach had substantial difficulties which indicate iiithat it did not have the flexibility to fully learn the data-set. Finally, the generic dense prediction network performed robustly in both the training and test data-sets. Additional testing is performed to contrast the performance of the two artificial neural networks.ivLay SummaryQuantifying data from micrograph images is often a laborious and time-consuming task involving manually tracing morphological features of cells. However, it is possible to train computational algorithms using manually traced data to automate image quantification, in a sense the machine learning algorithms learn by example how to perform the quantification task without additional human intervention. The goal of this thesis is to automate the quantification of myelin ultra-structure in transmission electron microscope images of the rat spinal cord. A variety of machine learning models are trained on images from spinal cord injury models producing myelin degeneration, attempting to quantify intact myelin area in images of different pathological and experimental conditions. The performance of the trained models is evaluated against manual tracing to determine their accuracy. It was found that a neural network based model performs this task in an accurate and unbiased manner, making it suitable to replace manual quantification of images. vPrefaceAll elements of this project involving animals were performed in accordance with the Canadian Council on Animal Care, and additionally approved by the UBC Animal Care Committee following Protocol #A12-0044. The initial project studying myelin debris using MRI and TEM, from which this paper derives its training image data, was suggested by Dr. Piotr Kozlowski. These training data have been used in a publication to address questions relating to how MRI scans relate to TEM ultra-structure (1). However, the analysis presented in this thesis is currently unpublished. The surgical injuries, animal care, sample preparation, and electron microscopy were done in Dr. Wolfram Tetzlaff’s laboratory at ICORD. The dorsal column transection surgical injuries were performed by Dr. Jie Liu. Animal care was performed by myself with the assistance of UBC staff veterinarians and ICORD-Jack Bell Centre animal care technicians. Tissue embedding, Histology, and imaging of TEM samples was performed by myself, and Dr. Henry Chen during a period when I was on medical leave (accounting for 35% of the TEM images used in this thesis). Manual tracing of images was performed by myself, Dr. Henry Chen, and two undergraduate volunteers under my supervision: Jason Wong, and Shawna Fu. TEM imaging was carried out under the supervision and guidance of Susan Shinn. All machine learning programming, training, and analysis was performed by myself.viTable of ContentsAbstract...................................................................................................................................................iiiLay Summary...........................................................................................................................................vPreface......................................................................................................................................................viTable of Contents...................................................................................................................................viiList of Tables............................................................................................................................................xiList of Figures.........................................................................................................................................xiiAcknowledgements.................................................................................................................................xvDedication..............................................................................................................................................xviChapter 1: Background Information....................................................................................................11.1 Cells of the central nervous system................................................................................................21.1.1 Neurons...................................................................................................................................21.1.1.1 Action potentials.............................................................................................................21.1.1.2 Electrical properties of neuronal membranes and cable theory of action potential conduction....................................................................................................................................41.1.2 Myelin and saltatory conduction............................................................................................51.1.2.1 Myelin structure..............................................................................................................61.2 The nervous system......................................................................................................................101.2.1 Central nervous system tissue composition..........................................................................101.2.2 Spinal cord............................................................................................................................11vii1.2.3 Dorsal column transection....................................................................................................121.3 Algorithms....................................................................................................................................151.3.1 Clustering.............................................................................................................................151.3.1.1 k-means, and Lloyd's algorithm....................................................................................151.3.1.2 k-Medoids.....................................................................................................................171.3.2 Classifying............................................................................................................................181.3.2.1 K-Nearest neighbour classification...............................................................................191.3.3 Neural Networks...................................................................................................................191.3.3.1 Artificial neurons and fully connected neural networks...............................................201.3.3.2 SegNet...........................................................................................................................251.3.3.3 Dense Prediction...........................................................................................................26Chapter 2: Machine learning image segmentation of intact myelin ultra-structure in TEM images rat spinal cord............................................................................................................................282.1 Introduction..................................................................................................................................292.2 Methods........................................................................................................................................342.2.1 Animal Preparation...............................................................................................................342.2.2 TEM fixation and imaging....................................................................................................362.2.3 General TEM image processing...........................................................................................382.2.4 General Image histogram processing....................................................................................382.2.5 Distance between image histograms.....................................................................................39viii2.2.6 Predicting sample condition from sample image histogram profile.....................................402.2.7 Semi-automated thresholding of myelin water spaces.........................................................412.2.8 Predicting intact myelin space in unsegmented images based on image histograms...........422.2.9 Computational Hardware......................................................................................................432.2.10 General neural network image preparation.........................................................................432.2.11 SegNet network training.....................................................................................................432.2.12 Dense prediction network training.....................................................................................452.3 Results and discussion..................................................................................................................472.3.1 Visualizing Histogram distances between and within data-sets...........................................472.3.2 Predicting sample condition from sample image histogram profile.....................................502.3.3 Semi-automated thresholding of myelin water spaces.........................................................632.3.4 Predicting intact myelin area in unsegmented images based on image histograms.............662.3.5 SegNet image segmentation and myelin area prediction......................................................712.3.6 Dense prediction image segmentation and intact myelin area prediction............................792.3.7 Expected error of machine learning approaches to quantifying sample myelin area...........852.3.8 Contrasting features of SegNet and dense prediction architectures.....................................87Chapter 3: Conclusion..........................................................................................................................913.1.1 Applications of non-parametric image histogram profile techniques...................................913.1.2 Applications of dense prediction for intact myelin segmentation........................................92Bibliography...........................................................................................................................................94ixAppendix.................................................................................................................................................99A.1 Complexity analysis..............................................................................................................99A.2 Big O, omega, and theta notation..........................................................................................99A.2.1 Practical assessment of time complexity..............................................................101A.3 Machine Learning................................................................................................................102A.4 k-means++ initialization......................................................................................................107A.5 Backpropagation of error in network training.....................................................................108A.6 Specialized network layers..................................................................................................111A.7 Image processing and Convolutional neural networks........................................................116A.8 Multi Dimensional Scaling..................................................................................................122xList of TablesTable 2.1 Confusion matrix of sample condition prediction for k1=8 k2=1, 3, 4, 5................................56Table 2.2 Linear model statistics for penalty parameters and their interaction.......................................62Table 2.3 Comparison of training and test correlation for prediction of semi-automated thresholding.. 65Table 2.4 Comparison of training and test correlation for prediction of unsegmented images by intensity histogram...................................................................................................................................70Table 2.5Comparison of training and test correlation for prediction of unsegmented images by SegNet. ....................................................................................................................................................75Table 2.6: Comparison of training and test correlation for prediction of unsegmented images by dense prediction.................................................................................................................................................82Table 2.7: Contrasting SegNet and dense prediction with modified architectures and training protocols. .................................................................................................................................................90xiList of FiguresFigure 1.1: Illustration of dorsal column transection injury.....................................................................14Figure 1.2: Representative TEM images of dorsal transection injury......................................................14Figure 1.3: Contrasting the components of a biological neuron and an artificial neuron........................24Figure 1.4: Sample artificial neural network...........................................................................................25Figure 1.5: Sample SegNet architecture..................................................................................................26Figure 1.6: Sample dense prediction architecture....................................................................................27Figure 2.1: Representative images of osmolarity experiment samples at different PB concentrations... 32Figure 2.2: Representative Image intensity histograms...........................................................................33Figure 2.3: Sectioning of spinal cord samples.........................................................................................36Figure 2.4: Within data-set MDS Sammon mapping for traced area histogram distances......................48Figure 2.5: Within data-set MDS Sammon mapping for whole image area histogram distances...........49Figure 2.6: Between data-set MDS Sammon mapping for traced area and whole image area histogram distances. .................................................................................................................................................50Figure 2.7: Exploring clustering and classifying parameters for predicting sample condition form imagehistogram profiles....................................................................................................................................55Figure 2.8: Sammon mapping visualization of histogram profile distances for k1=8.............................57Figure 2.9: Correlation between dispersion and accuracy for CVKNN classification of sample identity....................................................................................................................................................................58Figure 2.10: Accuracy of sample prediction with increasing CVKNN k2 values for selected k1 values..xii..................................................................................................................................................................59Figure 2.11: Accuracy of CVKNN sample prediction with increasing offset and size penalties............60Figure 2.12: Accuracy of CVKNN sample prediction with varying offset and size penalties for k2=1,3, and 5. k1=12 for each case.......................................................................................................................61Figure 2.13: Correlation and Bland-Altman plot of model training prediction vs. manual segmentation of intact myelin fraction for semi-automated thresholding......................................................................64Figure 2.14: Correlation and Bland-Altman plot of model test prediction vs. manual segmentation of intact myelin fraction for semi-automated thresholding..........................................................................65Figure 2.15: Correlation and Bland-Altman plot of model training prediction vs. manual segmentation of intact myelin fraction for unsegmented image prediction...................................................................68Figure 2.16: Histogram of injury data-set images by manually quantified myelin fraction....................68Figure 2.17: MDS Sammon mapping of untraced injury data-set image histogram distances coloured bymanual......................................................................................................................................................69Figure 2.18: Correlation and Bland-Altman plot of model test prediction vs. manual segmentation of image myelin fraction for unsegmented image prediction.......................................................................70Figure 2.19: Correlation and Bland-Altman plot of model training prediction vs. manual segmentation of image myelin fraction for SegNet myelin area segmentation..............................................................74Figure 2.20: Correlation and Bland-Altman plot of model test prediction vs. manual segmentation of image myelin fraction for SegNet myelin area segmentation..................................................................75Figure 2.21: Representative images of SegNet performance prediction of intact myelin in the training data-set.....................................................................................................................................................76xiiiFigure 2.22: Representative images of SegNet performance prediction of intact myelin in the test data-set.............................................................................................................................................................77Figure 2.23: Example image of SegNet prediction revealing errors in manual tracing...........................78Figure 2.24: Correlation and Bland-Altman plot of model training prediction vs. manual segmentation of intact myelin fraction for dense prediction myelin area segmentation................................................81Figure 2.25: Correlation and Bland-Altman plot of model test prediction vs. manual segmentation of image myelin fraction for dense prediction myelin area segmentation...................................................82Figure 2.26: Representative images of dense prediction performance prediction of intact myelin in the training data-set........................................................................................................................................83Figure 2.27: Representative images of dense prediction performance prediction of intact myelin in the test data-set...............................................................................................................................................84Figure 2.28: Expected myelin area sample error at increasing sample size for each machine learning approach...................................................................................................................................................86Figure A1: Graphical example of theta notation....................................................................................101Figure A3: Graphical example of Big-O and omega notation...............................................................101Figure A3: Demonstration of overfitting...............................................................................................106Figure A4: Sample network for illustrating error backpropagation.......................................................108Figure A5: Sample convolution of an image with a 2D filter................................................................112Figure A6: Sample max-pooling of an image........................................................................................113Figure A7 Sample unpooling layer.........................................................................................................114Figure A8: Example of an unpooling layer preserving low-level spacial information for recombination xivwith higher level processing...................................................................................................................115Figure A9: Convolution of an image with a horizontal filter.................................................................119Figure A10: Convolution as a constrained fully-connected network layer............................................120Figure A11: A simple convolutional neural network architecture for image classification...................121xvAcknowledgements I would like to express my gratitude to my supervisors, Dr. Wolfram Tetzlaff and Dr. Piotr Kozlowski, whose expertise, knowledge, and guidance I have greatly benefited from. Thank you for your patience with my medical problems, for giving me a 2nd, 3rd, and 4th chance, and for allowing me to take this thesis into a very odd direction.I would like to thank my additional committee member Vanessa Auld for her advice and feedback over the years at committee meetings (going back to my undergrad thesis).I would like to thank my collaborator Dr. Henry Chen for all his assistance, feedback, and mentor-ship during the initial experiments, and for keeping the project going when I was out of commission. I’m regretful that I left him high and dry as he was getting to the end of his thesis on this project.Special thanks to Susan Shinn for instructing me in the art and science of electron microscopy, and for keeping me company while I cut samples on the ultra-microtome in a literal closet for figurative weeks at a time.Thanks to Dr. Jie Liu for his expertise in surgery for the injury experiments, and for teaching me my way around a rat.Of course, I must thank my parents, and brother, who have supported me through all of this, that, and the other thing.Finally, I would like to acknowledge the Natural Sciences and Engineering Research Council of Canada and the Canadian Institutes of Health Research for funding this project.xviDedicationFor the strugglers.xviiChapter 1: Background InformationThe primary component of this thesis centres around the imaging and analysis of myelin in healthy tissue, and during degeneration. Several approaches to developing machine learning methods to automate the quantification of myelin area from transmission electron micrograph (TEM) images are assessed. This chapter will provide background information on the biology and structure of the myelin sheath, and an overview of the computational algorithms used in this work.The sections on cellular physiology and anatomy of the nervous system (§1.1-§1.2) provide background information on the function of myelin, and explains the anatomy of the injury model used to produce degenerating myelin in our experiment.The section on algorithms explains how the machine learning techniques used in this thesis function. The sections of clustering and classifying (§1.3.1-§1.3.2) describe the non-parametric approaches to analyzing electron micrographs based on their image intensity histograms which corresponds to the following methods (§2.2.6-§2.2.8) and results (§2.3.2-2.3.4) sections. The section on neural networks (§1.3.3) details the architecture of the networks used in this thesis (§1.3.3.2-§1.3.3.3) which correspond to the following methods (§2.2.10-§2.2.12) and results sections (§2.3.5-§2.3.6). Additional background information for readers unfamiliar with machine learning is contained in the appendix. §A.1-§A.2 explains some of the notation used to discuss computer algorithms. §A.3 Explainsthe basic idea of machine learning. §A.4 contains some additional details on initializing k-means. §A5-§A7 details the mechanics of artificial neural networks. Finally §A7 explains the algorithm used to visualize high dimensional data in this thesis.11.1 Cells of the central nervous system 1.1.1 NeuronsA neuron is the type of cell responsible for the processing and transmission of electrochemical impulseswithin the nervous system of animals. The word neuron derives from the Greek word for sinew or tendon, as both of these are tissues consisting of long thin fibres which have a similar appearance to neurons.A typical neuron consists of a cell body with two types of projections, dendrites and an axon. Many dendrites project a relatively short distance from the cell body; these projections receive signals from other neurons. Each neuron also has a single axon which may extend up to one meter in length, and may split at its end into many terminals allow the neuron to send signals to the dendrites of other neurons. The connection between the axon terminals of one neuron and the dendrites of the next are called synapses.1.1.1.1 Action potentialsA neuron at rest is kept polarized by sodium-potassium ATPase which moves three sodium ions out of the neuron in exchange for allowing two potassium ions to enter and is powered by the conversion of ATP to ADP. This results in a higher concentration of negative ions inside the neuron than outside, and the activity of many sodium-potassium ATPase pumps causes neurons to maintain a negative resting potential. This resting potential is used to allow neurons to conduct signals received from other neurons. When the synapses of a particular neuron’s dendrites are activated neurotransmitter-gated sodium or calcium ion channels open, causing the membrane of the dendrite to depolarize. Depolarization travels like a wave through the dendrites, and the cell body towards the start of the axon. Like a wave it gradually gets weaker as it spreads out, but if there is some combination of several2synapses activated at once, or one synapse activated repeatedly then the wave of depolarization will be strong enough at the start of the axon to trigger an action potential.An action potential is a strong and narrow wave of depolarization which is carried down the axon by the sequential opening and closing of voltage-gated ion channels along the length of the axon. A voltage-gated ion channel is a type of trans-membrane protein which changes conformation in responseto voltage changes across the membrane, opening a channel that allows certain ions to passively flow through it. An action potential is initiated when voltage-gated sodium channels at the axon hillock are depolarized past their activation threshold, which in turn causes them to open, allowing sodium to flow into the neuron along their concentration gradient causing further depolarization, which will in turn activate sodium gated ion channels further down the axon in a chain reaction. To allow the neuron to fire multiple action potentials rapidly the axon must be repolarized quickly after firing. To accomplish this the inactivation gate of the sodium-gated ion channels is closed which prevents further depolarization. The inactivation gate is triggered by the initial depolarization, but it takes some time to close so the sodium channel will be open for some time before it is blocked. The inactivation gate remains closed until the axon repolarizes. Additionally, the wave of depolarization also activates voltage gated potassium channels which allow potassium ions to flow out of the neuron along their concentration gradient causing the axon to repolarize. These potassium channels are activated at a higher membrane potential than the sodium channels, so temporally they will be activated after the sodium channels have been open for some time. A wave of sodium channels opening, being inactivated,and in turn activating potassium channels travels along the axon in step with the propagation of the wave of depolarization.When the action potential reaches the end of the axon it will cause the presynaptic-terminals to depolarize. This will trigger vesicles containing neurotransmitters to fuse with the membrane of the 3presynaptic terminal causing them to release their neurotransmitters into the synaptic cleft. These chemicals in turn bind to ion channels on the postsynaptic dendrite. This causes the ion channels to open and trigger depolarization, or repolarization in the postsynaptic neuron depending on whether the neurotransmitter released is excitatory or inhibitory respectively.1.1.1.2 Electrical properties of neuronal membranes and cable theory of action potential conductionThe speed and strength of action potentials as they travel down an axon depend on a number of factors which are best understood by abstracting the axon membrane as an electrical circuit. This involves dividing the axon into several subsections linked together in parallel, with each subsection having a resistor and a capacitor spanning the membrane. There are three properties that are important for this model. First the axon-length resistance rm which represents the resistance to current flow across the membrane for a short length of axon. This is derived from the area membrane resistance Rm, which is divided by the circumference of the axon Eq.(1.1)(2,3). Second the axon-length capacitance cm which represents the ability for the membrane to hold a charge across a short length of axon. This is derived from the area membrane capacitance Cm, which is multiplied by the circumference of the axon Eq.(1.2)(2,3). Finally the longitudinal resistance rl which represents the resistance to current flow within the cytoplasm of the axon along the length of the axon. This is derived from the specific resistance of the axon cytoplasm Rl, divided by the cross-sectional area of the axon Eq.(1.3)(2,3).rm=Rm2πa cm=Cm 2πa r l=Rlπa2 (1, 2, 3)These basic properties can be used to calculate two important properties of how the axon behaves in response to a change in applied voltage. The length constant λ represents how far an applied constant current can travel before the voltage change falls off to 37% of its initial voltage as it follows an 4exponential decay curve Eq.(1.4)(2,3). Increased membrane resistance will increase λ as this will ensure that more of the current flows along the membrane than across it, allowing it to travel along the axon more easily and increasing axon cytoplasm resistance will have the opposite effect. The time constant τ represents the amount of time it takes for an applied constant current to reach 63% of its steady state voltage after the voltage is applied following an exponential decay process Eq.(1.5)(2,3). Increasing both membrane capacitance and membrane resistance will increase the time it takes for the membrane to polarize. In the case of capacitance, more current is required to polarize the membrane, which makes the process of polarization take longer.λ=√ rmr l τ=rm cm (4, 5)1.1.2 Myelin and saltatory conductionRapid and efficient conduction of action potentials is important for animals to react quickly to their environment. There are two main adaptations which improve action potential velocity. One approach is increasing the axon's diameter. This decreases membrane resistance and increases its capacitance whichwould seem to make the axon's conduction properties worse Eq.(1.1, 1.2). However this also allows more paths for current to flow through the cytoplasm, which reduces longitudinal resistance by a greater margin Eq.(1.3). Overall increasing axon diameter increases action potential velocity linearly at a rate of ~5m/s/μm (4). The main drawback of this approach to improving action potential conduction is that axons will take up much more space, which can be limiting. Additionally thicker axons should require more energy for maintaining polarization, and in terms of production of proteins and other metabolic processes to maintain its basic structure.5The other approach to improving action potential conduction involves ensheathing the axon in a structure called myelin at regular intervals in order to isolate segments of the axon from the extracellular fluid. This helps prevent charge flowing between these spaces increasing membrane resistance, and spatially separates the charges of the spaces which reduces membrane capacitance. The effect of this is that myelination improves action potential velocity at a rate of ~15m/s for each micrometer of thickness added to the axon (5). This allows myelinated axons to remain narrower while still conducting action potentials quickly. The nature of action potentials travelling down a myelinated axon also changes. Instead of being sustained by a constant chain reaction of ion channel activity along the axon as in unmyelinated axons, the action potential quickly transmits through the myelinated sections of the axon to the unmyelinated nodes of Ranvier between segments of myelin. Ion channels are localized to these nodes along the axonand regenerate the action potential in the same manner as a normal action potential to allow it to speed through the next myelinated segment. This pattern of repeated potential regeneration and rapid travel along the axon is referred to as saltatory conduction. This pattern of action potential propagation also reduces the total amount of ion channel activity needed to conduct the signal, which in turn makes action potential conduction more energy efficient. 1.1.2.1 Myelin structureMyelin is produced by specific classes of glial cells within nervous tissue. In the peripheral nervous system axons are myelinated by Schwann cells which wrap themselves around the axon to produce myelin, requiring one Schwann cell per myelin sheath. In the central nervous system myelin is produced by oligodendrocyte which send out several projections of membranes to nearby axons, with each projection wrapping around the axon to produce a myelin sheath. This section describes the structure of oligodendrocyte myelin, general information on myelin structure is derived from the 6following sources unless otherwise cited (6,7).The structure of myelin is composed of a flattened membrane projection wrapped tightly several times around the axon. Were a myelin sheath to be unwrapped from an axon the area of compact myelin would roughly resemble an isosceles trapezoid with the projection from the oligodendrocyte attaching to the longer parallel edge, and the area contacting the axon being the shorter parallel edge. The ultra-structure of compact myelin is quite distinct in electron micrographs. Typical preparations stain the phospholipids and proteins associated with the membranes resulting in a darker stained major-dense band where opposed intracellular leaflets meet, and a fainter double-layered intra-period band where the opposed extracellular leaflets meet.The compact myelin is surrounded by a tube of loose membrane which is contiguous with the oligodendrocyte's cytoplasm. This space is referred to as the inner tongue, outer tongue, or lateral loop depending on whether it is located on the short parallel, long parallel, or angled edges of the trapezoid respectively. The general purpose of these spaces is to allow easier access for trafficking materials to compact myelin, and exchange of materials between the axon and oligodendrocyte. In cross-sections of compact myelin around an axon, myelin appears as a spiral of narrow extracellular and intracellular space separated by the membrane of the myelin projection. Immediately surrounding the axon is a small space of extracellular fluid and the inner tongue, while the outer loop is located outside of the compact myelin. At the distal ends of the compact myelin sheath are the paranodal regions where the lateral loops wrap around the axon in a helical pattern. The lateral loops are anchoredto the axon by structures called transverse bands. The combination of the repeated wrapping of compact myelin, and the anchoring of the lateral loops means that the extracellular space immediately surrounding the axon of compact myelin is effectively isolated from the rest of the extracellular space by a long narrow pathway.7The traverse bands are composed of tight junction proteins which hold the lateral loops firmly against the axon. For oligodendrocyte myelin these are composed of the proteins OAP-1/TSPAN-3 and OSP/claudin 11 (8). These proteins form a helical pattern in the paranodal region of myelin as they are aligned with the lateral loops which wrap around this region. At the outermost layer of myelin, the tightjunction continues along the outer tongue, connecting it to the outer layer of compact myelin in parallelto the ensheathed neuron, additional Tight junctions occur within the area of intact myelin parallel with the ensheathed axon, and aligned with each other radially between the layers of the sheath, known as the radial component of myelin (9,10). These structures help to restrict the diffusion of ions through theextracellular space of the myelin sheath improving its conductive properties (11).The compact myelin forms spontaneously as a result of intermolecular forces acting on it. Both the inner and outer leaflets of the myelin membrane contain lipids that feature exposed hydrophobic elements to favour driving water out, and hydrogen bonding elements to favour drawing opposed membranes closer together. Certain proteins are also trafficked into compact myelin. In oligodendrocytes the most common proteins are the transmembrane myelin proteolipid protein (PLP) and myelin basic protein (MBP). PLP functions based on the same intermolecular forces that membrane lipids use to stabilize myelin, however MBP has several basic groups which allow it to attract the negatively charged surfaces of the myelin membrane by electrostatic forces (12). The phospholipids and proteolipids of the myelin membranes also tend to have long and saturated fatty acidchains compared to other membranes (13). This makes compact myelin more densely packed and less fluid than other membranes. The combination of the effects of the structural features of myelin make it very stable. Knockout mice lacking individual myelin components often only exhibit pathological phenotypes in older age, and their myelin morphology often exhibits minor differences under normal conditions (14–16). The 8turnover rate of most myelin components are also extremely low, with half-lives beyond several months in radio-label experiments (17). The overall thermodynamic stability of myelin is such that intact compact myelin has been found preserved in 5000 year-old human remains (18). Note that the stability of compact myelin should not be taken to imply that it only functions as an inert material. Myelin also includes cytoplasmic pockets and gap junctions between adjacent layers of membrane (19).This allows transport between the axon and the oligodendrocyte through the layers of intact myelin, instead of requiring traversal of the spirals of the lateral loops. 91.2 The nervous systemIn mammals all sensory input, conduction, processing, and response is the responsibility of the nervous system which can be divided into two main components: the nerves and ganglia of the peripheral nervous system (PNS), and the brain and spinal cord of the central nervous system (CNS). The brain is encased by the skull and, in terms of communication with the rest of the nervous system, is responsible for processing sensory inputs and generating motor and endocrine responses. Impulses from the body are primarily conducted to and from the brain along the spinal cord. The spinal cord is a cylinder that extends in the spinal canal running from the brain to the upper lumbar levels in humans. It is responsible for some reflex arc and central pattern generation as well as integrating and conducting sensory and motor impulses between the brain and peripheral nervous system. The peripheral nervous system consists of bundles of nerve fibres which extend from the spinal cord at regular intervals, either as sensory afferent fibres from peripheral receptors and motor efferent fibres projecting to the muscles.1.2.1 Central nervous system tissue compositionWithin the CNS there are two main types of nervous tissue: grey matter which is responsible for processing information, and white matter which is responsible for rapidly conducting information. Bothtissues are distinctly separated spatially within the brain and are of significantly different composition. Grey matter is composed of neurons with a greater proportion of unmyelinated axons and are highly connected with other neurons in the vicinity in order to process sensory inputs. In contrast white matter is typically composed of parallel tracts of myelinated axons which allow quick conduction of motor, sensory and autonomic information over long distances. The myelin of the white matter is composed primarily of lipids-rich membranes wrapping around the axons of neurons in the tissue which gives it a somewhat lighter colour compared to grey matter.101.2.2 Spinal cordAs mentioned previously the spinal cord is one of the main routes for sensory/motor impulses between the brain and the peripheral nervous system. In terms of tissue the cylinder of the spinal cord is composed of an internal roughly butterfly-shaped grey matter (in cross sections) which is surrounded by white matter tracts. The grey matter of the spinal cord can be divided between dorsal, ventral, and lateral horns based on function. The dorsal horn of the spinal cord tends to be associated with sensory paths of information as sensory nerves enter the dorsal side of the spinal cord at regular intervals through the dorsal nerve rootsand some synapse onto the dorsal grey matter where the information can be processed to effect motor reflex output, or be sent to the brain through a white matter tract. Other sensory afferents (touch and vibration) enter through the dorsal root and send their processes directly through the dorsal columns to the brainstem nuclei where their 2nd order sensory neurons are located. The ventral horn of the spinal cord grey matter is in turn generally associated with motor pathways as it contains the cell bodies of motor neurons which project out of the spinal cord through the ventral spinal nerve root to reach the muscles which they innervate. Finally, the lateral horn mediates communication between the brain and visceral organs. neurons carrying descending autonomic impulses through the white matter hypothalamospinal tract synapse onto neurons of the lateral horn, which then conduct the impulse out through the ventral roots towards visceral organs. Similarly, sensory information from visceral organs enters the spinal cord through the dorsal root and synapse on neurons of the lateral horn. This circular connection allows some processing of autonomic impulses within the spinal cord. The grey matter is mirrored left to right such that the dorsal and ventral nerve roots for each side are responsible for sensor/motor impulses of their respective side. There is also a commissure of the grey matter across the centre of the spinal cord, allowing impulses to travel between the left and right sides 11of the spinal cord.The tracts of the white matter contain ascending sensory pathways and descending motor pathways. Of particular interest for this work is the dorsal-column of the spinal cord which contains three tracts in rats. The corticospinal tract descends from the motor cortex of the brain and carries motor impulses to motorneurons in the ventral horns at various levels of the spinal cord. In rodents it is the most ventral tract of the dorsal-column white matter and is of no interest in this study. The other two tracts of the rat dorsal column are the gracile fascicle and cuneate fascicle which carry sensory information from the lower and upper extremities respectively. The first order neuronal cell bodies of these fascicles are located in the dorsal root ganglion. These neurons have a branched axon that is split into two processes. The peripheral processes ends in sensory receptors that initiate action potentials, and the central processes enter the spinal cord through the dorsal horn to ascend through the spinal cord to the brain-stem where they synapse onto the 2nd order neurons in the gracile and cuneate nuclei respectively.1.2.3 Dorsal column transectionThe injury model used in my experiments involves a surgical transection of the dorsal column of the cervical spinal cord at the C5 vertebrate level. This particular injury model allows for a reproducible, discretely localized injury restricted to severing the white matter tracts of the dorsal column at a desiredlevel of the spinal cord. Caudal to the injury the axons remain attached to the neuronal cell body and survive the injury. However, the distal segments (cranial for the gracile fascicle) of the axon that are disconnected from their cell bodies begin to degenerate. Degeneration of the axons removes the survival signal for the 12oligodendrocytes in this area of tissue which causes their myelin sheaths to degenerate and form myelindebris that is gradually cleared by macrophages over a period of several months (20,21,21,22)Figure 1.1Figure 1.2. This is a process similar to Wallerian degeneration in the PNS where myelin begins disassembling structural components in a deliberate process. For example: endocytosis of tight junctionproteins, and fragmentation of compact myelin into droplets (23,24). Taken together this injury results in a consistent region of degenerating myelin sheaths that are spatiallyaligned with the gracile fascicle of the spinal cord cranial to the injury. This is a convenient model for studying the state of myelin during the process of degeneration. Additionally, a portion of the afferent neurons of the cuneate fascicle project caudally through the injury, which results in a combination of spared and degenerating myelin within the cuneate fascicle caudal to the injury (25). However this is ofno interest to the current study. Based on previous characterizations of this injury model, in which thereis no axonal regeneration, we expect to see only myelin debris in the severed tracts at 3-weeks post-injury, and partial clearance of these debris at 8-weeks post-injury (25)Figure 1.2.1314Figure 1.1: Illustration of dorsal column transection injury. Neurons of the gracile (red) and cuneate (green) fascicles enter the spinal cord and ascend through their respective tracts (solid lines). At the level of the transection injury the axons of the gracile are severed (grey region), The elements of these axons distal to the injury are no longerconnected to their cell body (dashed lines).Figure 1.2: Representative TEM images of dorsal transection injury. Images taken of the dorsal column 5mm caudal to the injury at 16 000x magnification. (left to right) healthy tissue, three weeks post-injury, eight weeks post-injury.1.3 AlgorithmsIn computer science an algorithm is a discrete and well-defined set of instructions which perform some specific task. Generally, an algorithm will act on some input in order to calculate a desired output. For example, a sorting algorithm will accept an array of numbers, process them through a series of steps, and finally return them in ascending or descending order. This section will present the general form of algorithms used within this thesis. Additional basic information is provided in the appendix which discusses some of the basics of algorithm analysis, general approaches to applying machine learning techniques, and mechanics of artificial neural networks.1.3.1 ClusteringThe purpose of a clustering algorithm is to divide a set of input points into several subsets such that the members of each subset are more related to one another than the members of other subsets by some metric. The typical input is a set of data-points, and the output is some metric which indicates which cluster each point belongs to. The specifics of the output will vary between different clustering algorithms as there is no single definition of a cluster. The specific algorithm used depends on the nature of the problem and the goals of clustering. 1.3.1.1 k-means, and Lloyd's algorithmk-means is an approach to clustering where the n input points are partitioned into k clusters Ck with each cluster being defined by a single point and each input being assigned to the cluster it is nearest to. An optimal clustering will minimize the within-cluster sum of squares distance Eq.(6)(29). Where nk is an element of n belonging to a particular cluster k, d is the number of dimensions defining each point orcluster, xnj is the value of variable j for point n, and xkj is the value of variable j for cluster definition k.15Error=∑k=1k∑n∈C k∑j=1d( xnj− x¯kj)2 (6)Practically, finding cluster definitions is usually approached using Lloyd's algorithm which accomplishes this by starting with seeded cluster centres, and iteratively refining them by finding all the points in n that are closer to ki than than all other cluster centres in k and updating ki with the averageof nk until the values of k fail to change after an iteration(30). Pseudo-code follows:While means in k have not convergedFor each na in nFind distance between na and each cluster mean in kAssign na to closest cluster mean in kFor each kb in kSet kb to be the average of all na assigned to kbThis approach to solving k-means does not guarantee finding the best possible cluster definitions, and the final result can depend on selecting good cluster seeds to avoid aberrant behaviour. These shortcomings can be mitigated by performing multiple runs and selecting the best result, and by using additional algorithms to find good cluster seeds. The advantage of this approach is that exhaustive search algorithms for solving k-means are computationally intractable as even simplified cases of k-means are in NP-hard (31,32). By contrast Lloyd's algorithm has much better performance in most cases. The Lloyd's algorithm approach to k-Means initially appears to have a time complexity of O(nkdi) where d is the dimensions of each point for distance calculation, and i is the number of iterations of the while loop before completion. This may appear to be effectively linear on n if the othervariables are fixed and/or much smaller than n for a given problem. However, the value of i is heavily dependent on the nature of the input data, sizes of n and k, and initial seeds. In the worst case Lloyd's 16algorithm will require 2√n iterations to converge (33), but for reasonable inputs it is likely to be of at worst some form of polynomial complexity(34).1.3.1.2 k-MedoidsA useful variant of k-means clustering is k-medoids. There are two substantial differences which distinguish k-Medoids. First, instead of cluster definitions being derived from averages of input points, only the input points themselves can be cluster definitions. Second, instead of attempting to minimize within-cluster sum of squares distance Eq.(6) K-Medoids selects cluster centres to minimize the sum ofpairwise distances between input points assigned to each cluster xn∈ Ck and the medoid of their clustersxk Eq.(7)(29). Where may be dist any desired distance function, k is the number of clusters. Error=∑k=1k∑n∈C kdist ( xn , x k) (7)There are several advantages to constraining the cluster centres to the input points in this approach. First, it avoids potential aberrant cases of k-means where the cluster definitions converge to spaces between clusters. Second, there are only n2 distances that can ever need to be calculated which allows them to be pre-computed. This allows any k-medoids algorithm to be entirely implemented in terms of matrix operations on the table of distances between each pair of points. Finally, K-medoids can be approached using a modified version of Lloyd's algorithm and k-means++, but the limited number of possible cluster centres makes greedy searches feasible and usually results in better clusters. A commonapproach to this is the Partitioning Around Medoids (PAM) algorithm (36). This algorithm starts with a set of k of the n data-points as medoids and attempts to swap each medoid with each non-medoid, maintaining swaps which reduce the sum of pairwise cluster distances Eq.(7), until the selection of medoids converges. Pseudo-code follows.17Determine quality of initial clusteringwhile medoids in k have not convergedfor each ka in kfor each nb in n, and not in kswap nb and kaassociate each n with nearest medoiddetermine quality of new clusteringif clustering is improvedupdate stored cluster qualityelseswap nb and kaThis algorithm has a time complexity of O(k(n-k)2i) where i is the number of iterations of the while loop before completion (initialization is of time complexity O(n2d), which typically will not be the dominant consideration). This may effectively give PAM a time complexity of O(n2) if k and i are muchsmaller than n. Similar to k-means clustering it is likely that this simplistic approach does not capture the influence of n, k, initial seeds, and the nature of the data on i. A paper published after the analysis inthis thesis was completed found that a modified version of k-means++ was effective at seeding k-medoids to improve clustering performance (37). Unfortunately a naive implementation of k-means++ was used to seed k-medoids in this thesis. 1.3.2 ClassifyingA classification algorithm is used to assort unknown points into one of several groups based on some metrics from a set of data. Similar to clustering there is a broad diversity of classifying algorithms which make different assumptions about the nature of the data, and are useful for different problems. Training a classifier requires having the training data-points separated into different categories, which 18is typically achieved through one of two approaches. The first is supervised learning where the category of each data-point is already known, for example when the data are derived from different treatment groups in an experiment. The second is unsupervised learning where data-points are not initially labelled, so categories must be discovered by first clustering the data. 1.3.2.1 K-Nearest neighbour classificationPerhaps the simplest classification algorithm is k-Nearest Neighbours (KNN) which works on the principle that unclassified points most likely belong to the same category as the known points that are closest to themselves. In its most basic form the algorithm works by counting the number of occurrences of each category in the k nearest points to an unknown point and assigning the category that is most common among those points (38). Classifying a single unknown point can be implemented in either O(n logn) time complexity for a sorting based approach, or O(nk) for a repeated minimum-find based approach. Where n is the number of known points. Clearly the minimum-find approach is faster when k<logn, but practically the difference is unlikely to be significant for most applications. One problem which complicates KNN is the occurrence of ties between multiple categories for highest frequency. This can be addressed by either ejecting the most distant neighbour in k until one category archives plurality, or weighting the contribution of each of the k neighbours by the inverse of their distance to the unknown point.1.3.3 Neural NetworksArtificial Neural Networks are an emerging area of machine learning research. The general approach involves building an architecture of several layers of artificial neurons, mimicking the abstracted structure and function of neurons in the human brain. These architectures typically consist of an input layer, multiple hidden layers where the inputs are processed by successive layers of neurons, and an 19output layer where the results of the processing are synthesized. Here we will explain the operation of abasic neural network, explain some types of layers that are used in this paper, show how collections of layers can be used to build a convolutional unit to allow image processing by convolutional neural networks, and describe the general forms of the neural network architectures used in this paper.1.3.3.1 Artificial neurons and fully connected neural networksAn artificial neuron can be imagined as a mathematical abstraction of a biological neuron. As discussedpreviously a neuron receives several synaptic inputs each of which may be inhibitory or excitatory to varying degrees based on synaptic plasticity (§1.1.1). If the sum of inputs is inhibitory no action potential will occur. If the sum of inputs is excitatory action potentials will occur with frequency proportional to the strength of the excitation. This functionality can be represented by extending our linear regression example from earlier to multi-variable linear regression, and passing the output through a rectified linear unit (ReLU).Multi-variable linear regression has a familiar loss function Eq.(1.8) and corresponding derivative Eq.(1.9) to single variable regression (27). However, x⃗ i and w⃗ now represent j-length vectors of sample features and weights respectively, where j is the number of variables. This makes a prediction yi for sample xi the dot product of multiplying each sample feature with it’s corresponding weight ( w⃗ x⃗ i = w1 xi1 +w2 xi2 +w3 xi3... wj xij), where xij is the jth feature of a sample i Eq.(1.10). This resembles the dendrites of a neuron where different synapses can have independently tuned excititory or inhibitory effects on the depolarization of that neuron. Similar to our earlier machine learning example §(A.3) gradient descent can be used with a data-set of dendritic inputs x, and corresponding depolarizing effecty to find a set of synaptic weights w which best explains how that neuron processes inputs.20f (w)=12∑i=1i(w⃗ x⃗i−yi)2 ∂ f (w)∂w j=∑i=1n(w⃗ x⃗i−y i) xij y i=∑j=1j(w j x ij) (8, 9,10)The ReLU in concept simply takes the maximum of zero and the result of the sum of weighted inputs to a neuron z Eq.(1.11) which represents the overall excitation level of the neuron(39). This models the conversion from neuron depolarization to rate of action potential firing. Either the sum of inputs prevents depolarization so that there is no action potential firing, or action potentials occur at a greater rate the more depolarized the neuron is. For the sake of using gradient descent for training a smooth approximation of ReLU is used in artificial neural networks Eq.(1.12) as it is convex and has a continuous derivative Eq.(1.13)(39). f (z)=max (0 , z) f (z)=ln(1+e z) ∂ f (z)∂ z= 11+e−z(11, 12, 13)This combination of a trainable multiple regression model and ReLU allows us to capture the general function of a biological neuron (Figure 1.3). There are some obvious abstractions made in this approach: synaptic tuning relies on both presynaptic and postsynaptic cellular machinery which is abstracted as a single postsynaptic regression weight in the artificial neuron. Discrete action potentials and time itself are also abstracted by expressing the activity of the neuron in terms of firing rate a fixed firing rate. Physical processes involved in basic neuronal physiology are non-linear; the length constantand time constant described earlier relating to membrane voltage changes describe exponential decay processes which can not be accurately captured by single linear weights§(1.1.1.2). Additionally dendritic signal integration has been found to have non-linear properties with regards to the processing of calcium action potentials, resembling the behaviour of an XOR logic gate(40,41). It is no coincidence then that representing the behaviour of a single simulated cortical neuron based on its observed inputs and resulting action potentials required a network consisting of five layers of artificial neurons as described here (42).21Now that we have our computational neuron assembled we can construct and train an artificial neural network. Typically, this will consist of an input layer of sample features, several hidden layers of neurons, and an output layer which contains one neuron corresponding to each possible sample classification. The layers are connected to one another such that the output of one layer becomes the input for the next layer: the input layer becomes the inputs for each neuron in the 1st hidden layer, the output of each neuron in the 1st hidden layer becomes the input of each neuron in the 2nd hidden layer, and so on until the output of the last hidden layer becomes the input of the output layer neurons (Figure1.4). In the output layer processing of neuron activation is performed differently to convert the set of neuron activation values to a set of classification probabilities for the corresponding sample identities by replacing the ReLU with a softmax layer Eq.(1.14)(27). Where Pi( x⃗ ) is the probability of a particular category out of n categories based on the vector of activations x⃗ from n neurons from the output layer of the network. This is an essential component for training as it constrains the output neurons such that their outputs range from zero to one, and sum to one. This can be interpreted as a probability for each possible classification. The gradient of an individual softmax probability Pi(x) withrespect to its direct activation activation xj where j=i, and other activations xj where j≠ i are respectivelygiven by Eq.(1.15, 1.16)(27).P i( x⃗)=ex i∑k=1ne x k ∂P i( x⃗)∂ x j=P i ( x⃗)(1−P i( x⃗))∂Pi≠ j( x⃗ )∂ x j=−P j( x⃗)Pi( x⃗) (14, 15, 16)The final prediction from the output layer is then compared with the expected sample identity to produce the error of prediction based on squared error Eq.(1.17) The gradient of the prediction error with respect to a particular category prediction is given by Eq.(1.18)(27).Error= 12∑i=1j(expected i−P i( x⃗))2 ∂Error∂Pi ( x⃗)=(P i( x⃗)−expected i) (17, 18)22There are several important considerations that arise in training this basic neural network using gradientdescent which are generally applicable to more complex variants. First the softmax layer is non-convexwhich means there is no longer a global minimum which gradient descent is guaranteed to reach, instead there may be many local minimums that the algorithm gets stuck in, and secondly the large number of weights in a neural network makes it computationally more expensive to perform an iteration of gradient descent. To address this a stochastic approach to gradient descent is used where foreach iteration or epoch the training samples are randomly divided into mini-batches which are then used to train the network separately. Introducing this randomness helps the network avoid getting stuck in local minima during training, and dividing the training data into mini-batches allows more adjustments to the model for the same computational effort. Training will typically be performed for many epochs while adjusting the learning rate (simply a multiplier to any adjustment to network weights) to gradually decrease how finely the network can be adjusted. When the learning rate is balanced properly the training process can be thought of as developing a general understanding for the trends in the data-set early in training, and having a precise understanding of uncommon cases, or complex relationships during late training. 2324Figure 1.3: Contrasting the components of a biological neuron and an artificial neuron1.3.3.2 SegNetSegNet is an approach to image segmentation using a convolutional neural-network architecture. The basic structure of a SegNet network is composed of an encoding half which uses 2D-convolution and pooling layers to develop a high-level understanding of the image, and a decoding half which uses max-unpooling layers and additional 2D-convolutions to apply that high-level understanding to segment the original image (50) (Figure 1.5). To conduct the unpooling in the decoding half of the network indices are used from the pooling half. This allows the pooling portion of the network to also encode boundary information at each level of resolution which can be combined with the higher level context during the decoding half by way of the trainable convolutional layers in the decoding half. The output layer has the dimensions of the original image with an extra dimension for the number of 25Figure 1.4: Sample artificial neural network categories for segmentation. The input to each pixel of the output layer is the filter response of that pixel from the preceding convolution (after batch normalization and rectification). During training eachpixel of the output therefore produces its own error which is then propagated backwards. 1.3.3.3 Dense PredictionThe dense prediction architecture used in this thesis is a different approach to image segmentation where conceptually a window is moved across the image, and the goal is to predict the central pixel of the window based only on the pixels within the window. For this study we are using a fully convolutional neural network with skip-connections that allow the output from shallower convolutionallayers to be passed forward to the final fully-connected layer which connects to the output layer (51) (To be clear SegNet is a specific kind of dense prediction network, I refer to the more generic architecture described here as dense prediction for convenience). The purpose of these skip-connections26Figure 1.5 Sample SegNet architecture. Network layers are scaled to show relative change in input dimensions. Lines connecting equal-dimension pooling and unpooling layers show the transmission of pooling indices. is to provide some low-level context to the final prediction and to make it easier for shallow levels of the network to have errors passed to them during training. 27Figure 1.6: Sample dense prediction architecture. Network layers are scaled to show relative change in input dimensions. The filtered image output of the first two pooling layers is fully connected to the output layers by a pair ofskip connections.Chapter 2: Machine learning image segmentation of intact myelin ultra-structure in TEM images rat spinal cord282.1 IntroductionThe original project that this thesis derives its training data from was intended to investigate a MRI technique for quantifying the presence of myelin sheaths based on T2 relaxation measurements. T2 relaxation can distinguish between water trapped in the layers of a myelin sheath and water in larger cell compartments, because water in the myelin sheath interacts with lipids in the sheath, causing it to relax faster than water in other compartments (53). This is used as a proxy measure to accurately quantify the fraction of myelin in healthy tissue, termed myelin water fraction (MWF) (25,53,54). In anearlier study involving MWF measurements performed on rats experiencing myelin degeneration following a dorsal column transection injury it was found that MWF measurements indicated an increased amount of myelin in the area of degeneration 3 weeks after injury (25,55). It was thought thatthis may indicate that as myelin degenerates more water is able to interact with the lipid membrane as itloses its structural integrity, which would account for the appearance of more myelin in the MWF MRI scan. To assess how ultra-structure changes of myelin sheaths correspond with MRI MWF measurements during myelin degeneration a series of experiments were performed comparing MWF results with manual traced vacuous, and intact myelin area in TEM images. It was found that the original anomalous MWF reading was the result of a fixation artifact during the processing of MRI samples, and that MWF otherwise correlates strongly with the quantity of manually traced intact myelin in TEM images (1,56). These experiments are the source of the training data-set in this thesis.One of the experiments performed while troubleshooting the above problem involved performing the MRI scans and TEM imaging in healthy samples at different concentrations of phosphate buffer to see if the changes in MWF or TEM morphology were sensitive to osmotic pressure, or variation in buffer concentration during fixation and MRI imaging (56,57). It was found that this intervention did not 29substantially affect the measured amount of intact myelin by MRI (after correcting for a shift in relaxation time), or traced TEM images. It did however, affect the morphology of intact myelin in the TEM images. Samples with buffer concentration substantially lower or higher than physiological conditions tended to have shrunken axons, and vacuous spaces in their intact myelin sheath, as if the sheath was no longer able to hold itself together, and started drifting apart Figure 2.1. This experiment was the source of the test data in this thesis.The set of TEM images, and corresponding manually traced myelin region masks produced for quantifying myelin ultra-structure in these experiments provides an opportunity to train a machine learning system to automate this quantification process. Machine learning techniques have recently been applied to different tasks in TEM myelin quantification (58–62). However, the higher level of magnification, and presence of myelin in pathological, or experimentally manipulated states in images produced for comparison between TEM intact myelin and MRI MWF necessitate training an approach to accomplish this specific task. In this thesis I attempt three separate approaches to automate the process of quantifying myelin in TEMimages. The first is based on KNN processing of image intensity histograms to infer the fraction of intact myelin in an untraced image based on the similarity of its intensity histogram to images with known intact myelin fractions. This is a similar concept to approaches for detecting sample misidentification in various data-sets (63,64), however it is typically used for classification rather than quantification tasks. The rationale for this approach is based on the observation that the intensity histogram for the traced myelin area is typically follows a bimodal distribution with the lower peak representing the stained myelin, and the higher intensity representing unstained water spaces (Figure 2.2 A, blue trace on graph). It is relatively straight froward to separate these peaks in the traced area only, but in the whole 30image there are other stained cell membranes, and unstained spaces. In the best case this will distort therelative heights of the traced area peaks peaks, and obscure the ideal intensity threshold to separate them (Figure 2.2 A, orange trace on graph). However, given that intact myelin stains darker than other membranes and proteins, these untraced membranes often fill in the space between the peaks in the whole image histogram, leaving a single skewed peak (Figure 2.2 B). Clearly there is too much noise relative to the signal we are interested in to measure the intact myelin in a whole image histogram directly. However, it may be the case that images with a similar intensity histogram have a similar intact myelin fraction. The other two approaches are two different artificial neural network architectures with contrasting priorities. The SegNet architecture attempts to learn a high-level interpretation of whole micrograph images that it is trained on; this serves as its basis for segmenting new images (50). By contrast the dense prediction architecture learns to segment images based on a low level understanding of how myelin appears in images, where only a limited region surrounding each pixel in an image is interpreted to segment that particular pixel (51). After training on the TEM injury data-set each approach will be tested on the osmolarity data-set to assess the performance of the model on an independently sampled data-set. The goal is to produce a model that can generalize the segmentation of myelin for use in future studies requiring the quantification of intact myelin area at high levels of magnification. It is important that a successful model is unbiased, and reasonably accurate relative to quantification by manual tracing.3132Figure 2.1: Representative images of osmolarity experiment samples at different PB concentrations. Images taken at 16 000x magnification. A) pilot sample images from samples fixed in (left to right) [0.0, 0..07, 0.14, 0.21, 0.28, 0..35] concentration of PB. B) samples from the experimental osmolatity samples fixed in (left to right) [0.023, 0.046, 0.07, 0.21] concentration of PB.33Figure 2.2: Representative Image intensity histograms. For each panel the top image is a TEM image, and the bottom image is the corresponding traced myelin area. The graph is the corresponding intensity histogram for the traced area (blue) and the whole image (orange). 2.2 Methods2.2.1 Animal PreparationAll experimental procedures were carried out in compliance with the guidelines of the CanadianCouncil for Animal Care and approved by the UBC Animal Care Committee prior to conductingthe study. All animals were housed in standard rodent care environments with food and water provided ad libitum. Of the nineteen male Sprague-Dawley rats that were used in this study twelve underwent the dorsal column transection surgery at the C5 level of the spinal cord as described by Chan et al (65). The rats were placed in a stereotaxic frame under anesthesia with their neck shaved and disinfected. Access to the spinal column was achieved by a midline incision along the axis of the spine followed by splitting the neck muscles at the midline to expose the cervical vertebrae. To access the spinal cord itself a laminectomy at the C5 level was performed using rongeur forceps. The transection of the dorsal column was carried out using a pair of micro-scissors pre-marked to indicate a depth of 1.8mm. A hypodermic needle was used to puncture holes on either side of the dorsal column to allow insetion of the scissor blades up to the marked depth, where the scissors were used to sever the dorsal column. Finally, to close the surgical site the neck muscles were sutured together followed by the skin.The injured animals were divided between 3-week and 8-week post-injury time-point groups (six animals in each), with six of the remaining animals comprising the uninjured group, and the final animal being used as a pilot for the osmolarity experiments. After surgery, the animals received an analgesic and were monitored frequently during recovery to ensure proper bladder function, nutrition, hydration, and excretion of waste. No significant injury complications occurred.34To sacrifice the animals, rats were deeply anaesthetized with isoflurane and perfused intracardially with0.1 M phosphate-buffered saline (PBS) for 3 min, followed by freshly hydrolyzed paraformaldehyde (4%) and glutaraldehyde (1%) in 0.1 M PBS at pH 7.4. The whole cervical spinal cord was then extracted and stored in a vial of the same fixation solution with the vials placed on ice for 30 minutes. For the main injury-based experiment spinal cords were sliced in the cross-sectional plane into a 1mm thick section for TEM and a 3 mm section for MRI, with the TEM section cranial to the MRI section and the division between the two sections being centred 5 mm cranial to the injury site (Figure 2.3). The 3 mm MRI sections were further fixed post-extraction overnight in 2% glutaraldehyde before scanning, while the 1 mm TEM segments were fixed immediately TEM imaging (next section). These 1 mm TEM segments collectively make up the training samples for the present study.The non-injured rats had additional segments removed for a separate osmolarity-based experiment. For each animal the remaining cervical spinal cord was sliced into six additional sections, alternating between 1 mm and 3 mm sections with each pair of adjacent 1 mm and 3 mm sections being a matched pair for TEM and MR imaging respectively (Figure 2.3). Each pair of samples was fixed in a 2% glutaraldehyde solution of one of a set of PB concentrations: 0.21 M, 0.07 M, 0.046 M, 0.023 M, and 0 M. Samples were randomly assigned with the constraint that no one animal had more than one sample pair fixed at a given concentration, and samples of a particular PB concentration were taken from different levels of the cord where possible. The osmolarity pilot animal was treated similarly, however it had six sample pairs removed and fixed at the following PB concentrations: 0.35 M, 0.28 M, 0.21 M, 0.14 M, 0.07 M, and 0 M. All sample pairs were fixed overnight before being separated for MR and TEM fixation (next section). These 1mm TEM segments collectively make up the test samples for the present study. 352.2.2 TEM fixation and imagingFor the duration of the fixation processing the 1mm spinal cord samples were kept on ice to maintain temperatures close to 4°C. Sections were postfixed for one hour in 2% glutaraldehyde for lipid and protein fixation. The sections were then washed in 0.2M cacodylate buffer for 3×10 minutes and placedin 1% osmium tetroxide and 1.5% potassium ferrocyanide in 0.1M cacodylate buffer for one hour for lipid fixation. After fixation samples were rinsed for 10 minutes in Milli-Q purified water then dehydrated through rinsing with a series of increasing concentrations of EtOH (2×5 minutes in 50% EtOH, 2×5 minutes in 70% EtOH, 3×10 minutes in 95% EtOH, and 3×10 minutes in 100% EtOH). Sections were rinsed in propylene oxide (PO) for 3×20 minutes then Spurr resin was infiltrated into the tissue through a series of rinses (1 PO : 1 Spurr overnight, 1 PO : 3 Spurr 3 hours, 100% Spurr 6 hours, 36Figure 2.3 Sectioning of spinal cord samples. Injury site (red), injury experiment sample pairs (orange), and osmolarity experiment sample pairs (green) are shown relative to the C5 level of the spinal cord.100% Spurr overnight). Sections were then embedded into blocks using fresh 100% Spurr resin which was polymerized for 24 hours at 60°C.Once embedded the samples were sectioned on a Leica Ultracut T ultra-microtome (Leica Microsystems AG, Wetzlar, Germany) using 45° glass knives. 1μm thick semi-thin sections were made of the entire cross-section of each sample to determine the area of the dorsal column in each sample before cutting the blocks down to that area, and to serve as a guide during TEM imaging. These sections were stained with 0.5% toluidine blue and were imaged with a Zeiss Axioplan 2 microscope (Carl Zeiss AG, Oberkochen, Germany) at 200× magnification. Once the blocks were trimmed to the area of the dorsal column 70nm ultra-thin sections were cut using the same microtome setup as above. These sections were mounted on copper SPI TEM grids (85×85μm windows, 40μm thick grid bars) then counterstained with 2% uranyl acetate for 10 minutes at 60°C, followed by Reynolds’ lead citrate for 2 minutes in the presence of sodium hydroxide pellets. All TEM images were taken on a Zeiss EM910 transmission electron microscope (Carl Zeiss, Thornwood, NY) using an AMT digital camera (American Technologies Corporation, San Francisco, CA) and iTEM analysis software (Olympus Soft Imaging Solutions GmbH, Münster, Germany, version5.0). For each animal three grid windows in the gracile fascicle were located for sampling. Nine images were taken at 16,000× magnification, in a 3×3 fixed pattern such that ~10μm separated each 5.83µm × 5.83µm image. Images were captured and stored at 2048×2048 pixels resolution, with intensity levels quantized to 16-bit levels. One animal from each of the 8-week post-injury and 3-week post-injury groups were excluded from analysis due to poor sample fixation. Additionally six of the 18 osmolarity experimental samples had poor fixation resulting in the following number of samples at each phosphate buffer concentration: one at 0.21 M, two at 0.07 M, four at 0.046 M, three at 0.023 M, and two at 0 M.372.2.3 General TEM image processingThe intensity values for the 16,000× magnification TEM images were rescaled by mapping the 0.1 to 99.9 percentile pixel intensities of each image into the 0 to 1 range linearly, with resulting values belowzero or above one being set to zero or one respectively. To reduce random noise and balance uneven image illumination an 8×8 adaptive Wiener filter, followed by contrast-limited adaptive histogram equalization (8 × 8 tiles and 256 bins) were applied to the images. This resulted in 2048×2048 pixels resolution images with 8-bit depth intensity levels. The rightmost 136 pixels of each image were excluded from analysis due to a persistent TEM imaging artifact in the form of a shadow over the right edge of several samples, giving a final sample resolution of 2048×1912, To isolate the area of myelin sheaths each image was manually segmented by tracing the area of both intact and degenerating myelinsheaths. Each segmented image was saved as a separate binary image. To find the optimal thresholding for separating intact myelin from vacuous myelin spaces in the segmented myelin space each image was initially thresholded using Otsu's method (66). This initial thresholding was manually inspected using an interactive image thresholded graphical user interface (GUI) in MATLAB, and manually adjusted in cases where the threshold did not correctly contrast intact myelin from vacuous spaces.2.2.4 General Image histogram processingThe image histograms of the TEM images captured from all samples were processed to extract the intensity histograms for the histogram-based experiments. For both the myelin area and the whole image area the raw histogram was smoothed with a one dimensional Gaussian filter (σ = 2 and size =7).The resulting intensity histograms were quantized to 100 intensity levels, and normalized by dividing by the myelin or image area respectively. 382.2.5 Distance between image histograms.For all histogram experiments the distance between two intensity histograms was calculated based on the following equations:FinalDist= min−25≤i≤25SizeCost+(|i |×OffsetPenalty /2500)+∑j=1100(histA j−histB j+i)2100−|i|(19)Where histAj is the value of the jth point of a particular histogram. The rightmost summation represents taking the sum of the square distance between each pair of points on the histograms. Since the goal is tocompare histograms based on their general shape we allow the histograms to be offset by up to 25 bins in each direction using the variable i. The denominator normalizes the distance of an alignment based on the number of overlapping points in that alignment. As we want to compare histograms at their best alignment we will accept the minimum distance alignment as the FinalDistance. The OffsetPenalty is a value greater than, or equal to, zero. This is multiplied by the magnitude of i to apply a scaling penalty to each alignment. SizeCost is an additional penalty based on the difference in myelin area traced whichis calculated based on the following equation:SizeCost=[ |log10(1+myea/ tota)−log10(1+myeb/ totb)|×SizePenalty /100] (20)Where myea is the segmented area of myelin in a particular image, tota is the total quantifiable area of that image, and SizePenalty is a value greater than, or equal to, zero. This applies a penalty based on the difference in relative myelin area of the image.The arbitrary constants in these formulas were a practical work-around to move the dynamic range of penalty values to be whole integers. Through experiments in this section the OffsetPenalty is included as a hyper-parameter given the possibility that highly offset alignments should only be accepted if they strongly improve the distance 39score. This may help avoid aberrant alignments. The SizePenalty variable is included as a hyper-parameter given that intensity histograms of images with vastly different traced myelin areas are possibly different in their characteristics in ways that are not fully captured in the intensity histogram (essentially the size penalty partially ‘un-normalizes’ the histograms for the traced area histogram experiments). To put these parameters in context the average pairwise histogram distance for all images in both data-sets is 0.0470. An OffsetPenalty of 10 therefore adds a penalty of 0.004, or 8.5% per offset on average. Similarly, the average pairwise distance for the log component of SizeCost is 0.601. Therefore a SizePenalty of 10 adds a penalty of 0.0616 or 106% on average.This approach for calculating distance between histograms introduced some constraints on the machinelearning and algorithms used for analysis of this data. As histograms are compared with different offsets their distances do not represent distances between points in euclidean space (distances between three histograms may violate triangle inequality, such that the longest euclidean distance between the three histograms is longer than the sum of the two shorter distances). This means that the k-means update step is attempting to find an average point within non-euclidean space which is not guaranteed to reduce error of clustering, and often failed to converge. k-medoids improves on this as it does not attempt to define an average point in non-euclidean space, and PAM ensures that the error of clustering will decrease monotonically (36).2.2.6 Predicting sample condition from sample image histogram profileThe TEM images from the injury-based experiment were clustered based on the distance between image histograms of their manually traced area using k-medoids clustering with a euclidean distance function. This was repeated with a range of k1-values from one to twenty, with thirty repetitions per value to attempt to find a value with the best average clustering performance.40For each individual sample an image histogram profile was found by determining the proportion of images from that sample that occur in each cluster as determined by k-medoids clustering with a euclidean distance function. This produces an array of k1 values between zero and one for each sample. Each sample's identity was predicted by comparing image histogram profiles among all samples using k-fold cross-validated k-nearest neighbour classification (CVKNN). The range of k2-values were attempted from one to ten with thirty repetitions per value to maximize prediction accuracy.With suitable k1 and k2 values found, clustering and classifying was repeated using the selected k-valuesto optimize the hyper-parameters OffsetPenalty and SizePenalty. OffsetPenalty was varied from zero to twenty-seven by steps of three and SizePenalty was varied from zero to 100 by steps of ten. Sixty clusterings were attempted per combination, with each sample's identity predicted CVKNN to measure validation accuracy.2.2.7 Semi-automated thresholding of myelin water spacesThe TEM images from all sample conditions of the injury-based experiments were clustered based on the distance between manually-traced myelin area intensity histograms using k-medoids clustering. To predict the optimal thresholding for segmenting intact myelin from vacuous myelin spaces the manually determined threshold point of each medoid was mapped onto the corresponding point of eachmyelin intensity histogram in its cluster, while taking into account how each histogram alignment was offset from the medoid. The resulting intact myelin area after thresholding was divided by the area of the image to obtain the myelin fraction for each non-medoid image. To optimize the model’s predictive accuracy the above process was repeated at several different k1, OffsetPenalty and SizePenalty values using greedy descent to select ideal parameters to minimize total squared prediction error of myelin fraction compared to manual tracing and thresholding. Error for 41each image was found using k-fold cross-validation. The k1 value was varied from one to twenty, OffsetPenalty was varied from zero to twenty-seven by steps of three, and SizePenalty was varied from zero to 100 by steps of ten.With an appropriate k1 parameter, and the hyper-parameters OffsetPenalty and SizePenalty found the final model was tested by predicting using the same parameters to perform the semi-automated thresholding process to determine the myelin fraction of the TEM images from the osmolarity samples.2.2.8 Predicting intact myelin space in unsegmented images based on image histogramsThe TEM images from all sample conditions of the injury based experiments had their unsegmented and untraced intensity histograms and manually determined myelin fraction data combined to produce a model to predict the myelin area of new images. Prediction is accomplished by a variation on KNN where the myelin fraction was estimated by taking the mean myelin fraction from the k-nearest image histograms in the model. Validation accuracy of the model was measured by k-fold cross-validation to compare the manually measured myelin fraction with the predicted myelin fraction for each sample in turn. To optimize the model’s predictive accuracy the above process was repeated at several different k2, and OffsetPenalty values using greedy descent to select ideal parameters to minimize total squared prediction error of myelin fraction compared to manual tracing and thresholding. Error for each image was found using k-fold cross-validation. The k2 value was varied from one to twenty, and OffsetPenalty was varied from zero to twenty-seven by steps of three.With an appropriate k2 parameter, and the hyper-parameter for OffsetPenalty found, the final model wastested by predicting the myelin fraction of the TEM images from the osmolarity samples.422.2.9 Computational HardwareAll training was performed on a home computer primarily composed of an Intel i5-6600K processor (Intel Corporation, Santa Clara, CA), NVIDIA GeForce GTX 1060 6Gb VRAM video card (Nvidia Corporation, Santa Clara, CA), and 16Gb of RAM.2.2.10 General neural network image preparationAll images were shrank by a factor of eight to a resolution of 256×239 for the neural network experiments (effectively resembling the resolution of a 2000x magnification image). This was necessary to fit the training data-sets in memory and helped reduce training time by reducing the size ofnetworks. 2.2.11 SegNet network trainingA five-level deep SegNet architecture was prepared for this experiment. Each of the five sequential levels on the encoding half of the network was composed as follows:• 2D-Convolutional layer: 3x3 field, 64 filters• Batch normalization layer• Rectified linear unit• 2D-Convolutional layer: 3x3 field, 64 filters• Batch normalization layer• Rectified linear unit• Max Pooling: 2x2 field, 2,2 strideWhile the decoding half was similarly composed of five sequential levels each as follows:• Max-Unpooling layer43• 2D-Convolutional layer: 3x3 field, 64 filters• Batch normalization layer• Rectified linear unit• 2D-Convolutional layer: 3x3 field, 64 filters (2 for final output)• Batch normalization layer• Rectified linear unitWhere each max-unpooling layer is taking indexes from the max pooling layer of corresponding level in the encoding half. The final rectified linear unit feeds into a softmax layer to produce the network output. The input to the network is the 256×239 TEM micrograph and the corresponding training targetis the binary image of manually segmented and thresholded myelin area. The TEM images from all sample conditions of the injury based experiments were used to train the network. The training protocol consisted of 300 epochs with a minibatch size of fifteen. Training samples were shuffled and augmented every epoch by randomly rotating them by up to 360° and offsetting them by ±10 pixels vertically or horizontally, while performing the same transformation to the output image. Training was conducted at an initial learning rate of 0.01 which was dropped every 40 epochs by a 0.1 multiplier. All weights were regularized by the square of their magnitude multiplied by a 0.0001 regularization weight.During training Validation error was checked every 300 mini-batches by measuring the accuracy of the network in predicting the original non-augmented images of the training data. This was done rather than splitting the training data into a separate validation data-set due to the limited number of images inthe training set.44After training the network was used to predict the manually segmented and thresholded myelin images for the osmolarity experiment TEM images. 2.2.12 Dense prediction network trainingA five-level deep dense-prediction architecture was prepared for this experiment. Each of the five sequential levels consisted of a convolutional unit as follows:• 2D-Convolutional layer: 3x3 field, 32*2i filters• Batch normalization layer• Rectified linear unit• Max Pooling: 2x2 field, 2,2 strideWhere i is the level number. The 1st,3rd, and 5th layers had their outputs fed into a dropout layer then fully connected to a softmax layer, which output the two categories. The category prediction from each of the three fully connected layers was combined by an additional layer which was fed into a softmax layer to produce the output. The input to the network is a 32x32 pixel patch of a TEM micrograph and the corresponding training target is the centre pixel of the patch from the binary image of manually segmented and thresholded myelin area. The TEM images from all sample conditions of the injury based experiments were used to train the network. 399 non-overlapping patches were extracted from each of 523 images. Each patch was also rotated by 0°, 90°, 180°, and 270° to augment the training data-set. This resulted in a total of 834 708 patches; 10 000 were randomly selected to use as a validation data-set during training. Training was performed for 35 epochs with a minibatch size of 500. The initial learning rate was 0.003 and was dropped every 6 epochs with a drop factor of 0.2. Validation accuracy was assessed every 500 iterations during training. All weights were regularized by the square of their magnitude multiplied by 45a 0.0001 regularization weight with the exception of the fully connected layers that came from the 1st and 3rd levels of the network where the regularization weight was increased to 0.001. Dropout probability was set to 50%.After training the network was used to predict the manually segmented and thresholded myelin images from the osmolarity experiment TEMs by applying the network to predict each pixel one-by-one.462.3 Results and discussion2.3.1 Visualizing Histogram distances between and within data-setsMulti-dimensional scaling visualizations of histogram distances for both the injury and osmolarity data-sets were produced using Sammon mapping to aid in intuition (29). Histogram distance matrices were calculated with both the SizePenalty and OffsetPenalty parameters set to zero for the visualization. The within data-set visualizations for the traced image area histograms show some segregation of the distribution of histograms by sample categories in the injury data-set, and gradual shift in the sample category distribution with increasing buffer concentration for the osmolarity data-set (Figure 2.4). Similar trends exist in the whole-image histogram distance visualization where the shift in the osmolarity data with respect to concentration is more pronounced (Figure 2.5). These visualizations give an indication that the image intensity histograms of sample images are influenced by sample condition both within the myelin space, and for the whole image area. This supports the histogram-distance based machine learning approaches applied in this section. A visualization of between data-set histogram distances was produced for both traced area and whole image histograms (Figure 2.6). Histogram samples from the osmolarity data-set appear within the area injury data-set, with some differences in density between the two data-sets. This gives an indication that a machine learning model trained on the injury data could be applicable to the osmolarity data-set as samples from the two data-sets are not segregated. 4748Figure 2.4: Within data-set MDS Sammon mapping for traced area histogram distances. Injury (top) and osmolarity(bottom) data-sets are visualized using Sammon mapping. Points from each data-set is plotted based on their sample condition. Axes are unitless dimensions.49Figure 2.5: Within data-set MDS Sammon mapping for whole image area histogram distances. Injury (top) and osmolarity (bottom) data-sets are visualized using Sammon mapping. Points from each data-set is plotted based on their sample condition. Axes are unitless dimensions.2.3.2 Predicting sample condition from sample image histogram profileA manual approach to parameter selection and model optimization will be used for predicting sample conditions to examine the performance of histogram-distance clustering based approaches used in later 50Figure 2.6: Between data-set MDS Sammon mapping for traced area and whole image area histogram distances. Traced area (top) and whole image area (bottom) distances are visualized using Sammon mapping. Points are plotted based on their data-set. Axes are unitless dimensions.sections for myelin area prediction models. Average dispersion was measured for thirty repeated clusterings of the injury data-set image histogramswith k1-values for medoid clustering set from 1 to 20 with a euclidean distance function. Approaches to identifying a suitable k1 statistically using the gap statistic were unsuitable given that the points of the intensity histogram are non-independent, and there is no obvious way to map the distances between image histograms into Euclidean space (67). By inspection the dispersion of clustering attempts ceased to improve substantially for k1 ≥ 6 (Figure 2.7 A). This suggests that there are six different 'types' of image histograms in the data. All clusterings were carried forward to perform an initial attempt at CVKNN to predict sample conditions with k2 = 1, 3, 4, and 5. The accuracy of classification improved from k1=1 to k1=5 after which changes in accuracy were less drastic (Figure 2.7 B). The k2=1,3,4 cases continued to improve slightly until k1=8 after which they levelled off, while k2=5 had slightly worse accuracy in comparison. The lack of accuracy when k1<5 is likely due to some of the useful signal for predicting sample conditions based on the profile of image histograms being combined into the same clusters, which as a result causes histogram clusters to not segregate between sample conditions. Withinthe repetitions of a single k1 and k2 value with good accuracy there was no correlation between the dispersion of clustering, and the accuracy of classification, so classification performance can not be improved by repeating clustering several times and performing CVKNN on the least dispersion clustering (Figure 2.9).The confusion matrices of the k1=8 k2=1, 3, 4, and 5 sets of classification repetitions indicates that classification is fairly robust (Table 1). The largest source of the error in prediction is the result of a single uninjured sample being consistently miss-classified as an 8-week post-injury sample, similarly there is a 3-week injury sample often miss-classified as an 8-week sample, and an 8-week sample that is often mistaken as uninjured. This suggests that this particular sample is a borderline case that is 51difficult to classify given the nature of the data, rather than the clustering and classification processes themselves being generally inaccurate. A Sammon mapping visualization based on histogram profile distances for k1=8 shows that samples generally cluster based on sample identity, and that the particularuninjured sample that is frequently sampled (labelled Nor3) is in the space between the uninjured and 8-week post-injury sample clusters (Figure 2.8). Different clusterings of the image histograms may move this sample’s histogram profile closer to either cluster which explains why it is consistently miss-classified. Similarly, the 8-week post-injury sample that is often miss-classified as uninjured (8wk1) is separated from the main cluster of 8wk samples and is actually closer to the Nor3 sample. It is less clear from this visualization why the 8-week post-injury sample is often misclassified as a 3-week sample (3wk4), certainly it is somewhat separated from the other member of its sample group, and is closer to the border between sample categories, but it is also much closer to members of its own samplecategory than the 8-week post-injury cluster. We know it should be nearly equidistant based on the error rate, clearly there is some distortion when we are mapping points from an 8-dimensional space into a 2D figure.The k2=1 case is of particular interest as a lower bound for the Bayes error of classification of the data-set (68). The maximum accuracy values for k2=1 is occasionally 100% so the current data-set can be clustered such that the image histogram profiles perfectly segregate based on sample identity. However,for larger sets of samples one expect that the underlying distribution of image histogram profiles wouldresult in a non-zero Bayes error of classification, given that several samples already have a tendency to be incorrectly classified (Table 1).Going forward a k1=8 and 12 were selected to be used for image histogram clustering by inspection, as k1=8 seems to be where performance improvement begins to level off, and k1=12 is the lowest value beyond which there is no apparent improvement in accuracy. To determine suitable k2 values for 52classification CVKNN was repeated with k2 values from 1 to 10, for 60 image histogram clusters each (Figure 2.10). The accuracy of the classification was greatest at lower k2 values, and relatively stable for k2≤4 for both k1 values. At greater k2 values accuracy declined sharply. There were only 5-6 samples for each treatment condition, so performing CVKNN with increasing k2 guarantees that an increasing number of incorrect samples will influence the classification of any particular point.To examine the effect of Offset penalty, or Size penalty adjustments to the distances between image histograms on classification and clustering, K-medoids and CVKNN was repeated with k1=12, k2=1, 3 and 5, and both penalty values being stepped independently from zero to 100 by steps of 10 for Size penalty, and from zero to 30 by three for Offset penalty. These dynamic ranges were identified by manual inspection and adjustment starting from a max penalty of 1000, and working downwards. When varying each penalty independently Size penalty tended to improve sample prediction accuracy significantly and with a positive trend, for all k2 values (Figure 2.11). This does not directly indicate that the amount of traced area in an image depends on sample condition, but that applying this penalty to histogram distances produces a k-medoids clustering that results in image histogram profiles that are in turn better clustered by sample category. In contrast, Offset penalty tended to reduce sample classification accuracy significantly and with a negative trend for k2 = 3,5 while improving accuracy significantly for k2 = 1 for Offset penalty = 3,6 and having no effect for larger value (Figure 2.11). Given that we know the k2 = 1 case is representative of the Bayes error rate we can infer that increasing Offset penalty still leaves the sample category clusters for the histogram profiles distinct from one another as we saw earlier in the MDS visualization. However, given that the k2 = 3,5 accuracy decreasessubstantially we can infer that the distances between samples within a sample category cluster is increasing relative to the distances between sample category clusters. 53To see if how the combination of the Offset penalty, and Size penalty together affects accuracy, both were varied together for the k2=1,3, and 5 data to generate a heat map (Figure 2.12,Table 2). The top row and left column are the same data as the previous figure (Figure 2.11). The interaction of these parameters appears independent of one another for k2=3,5: increasing Offset penalty reduces accuracy, and increasing Size penalty improves accuracy. This is reflected in the linear model fit to these data as the interaction variable is positive, but has an estimated coefficient two orders of magnitude smaller than the coefficients of the parameters themselves, and the interaction coefficient for k2=3 case does not meet the p<.05 level of significance, suggesting that it is not distinct from a coefficient of zero. For the k2=1 case the Offset penalty coefficient itself is not statistically distinct from zero, which is congruent with varying the Offset penalty not impacting classification accuracy as we saw previously. The interaction coefficient here is also small but negative. Overall this linear model shows that the penalty parameters influence the classification accuracy largely independently from one another across k2 values tested.Altogether in this section we have seen that image intensity histograms have some relation to a concrete property of our samples, as the image histogram profiles for our samples are able to predict sample identity. Additionally, we have explored the manual optimization of parameters for k-medoids clustering, k-nearest neighbours classification, and the particular penalty parameters for the image histogram distance calculation. 5455Figure 2.7: Exploring clustering and classifying parameters for predicting sample condition form image histogram profiles. A: Average dispersion resulting from the clustering of the set of images at by k-medoids at increasing k1 values. B: Average accuracy of sample condition prediction by CVKNN for various k2 values at increasing k1 values for the initial clustering (solid lines with dots). The best accuracy of the k2 =1 attempts are also plotted (dashed lines). Error bars are standard error of the mean in all cases.56Table 1 Confusion matrix of sample condition prediction for k1=8 k2=1, 3, 4, 5. Average frequency of CVKNN condition prediction for each sample across 30 repetitions are shown. Green cells represent the expected sample identity. Orange cells represent samples that were incorrectly classified in more than half of repetitions, yellow represents incorrect classification in over a quarter of repetitions.Classifier Sample IdentityParameter Nor1 Nor2 Nor3 Nor4 Nor5 Nor6 3wk1 3wk2 3wk3 3wk4 3wk6 8wk1 8wk3 8wk4 8wk5 8wk6Predicted IdentityNor 1.00 1.00 0.47 1.00 1.00 0.93 0.00 0.00 0.00 0.00 0.02 0.52 0.00 0.00 0.00 0.003wk 0.00 0.00 0.07 0.00 0.00 0.02 1.00 0.92 0.90 0.58 0.98 0.08 0.15 0.00 0.02 0.088wk 0.00 0.00 0.47 0.00 0.00 0.05 0.00 0.08 0.10 0.42 0.00 0.40 0.85 1.00 0.98 0.92Nor 1.00 1.00 0.40 1.00 1.00 0.88 0.00 0.00 0.00 0.00 0.00 0.20 0.00 0.00 0.00 0.003wk 0.00 0.00 0.05 0.00 0.00 0.05 0.93 0.82 0.92 0.57 1.00 0.08 0.17 0.00 0.00 0.028wk 0.00 0.00 0.55 0.00 0.00 0.07 0.07 0.18 0.08 0.43 0.00 0.72 0.83 1.00 1.00 0.98Nor 1.00 1.00 0.40 1.00 1.00 0.90 0.00 0.00 0.00 0.00 0.00 0.22 0.00 0.00 0.00 0.003wk 0.00 0.00 0.05 0.00 0.00 0.05 0.93 0.82 0.92 0.58 1.00 0.10 0.17 0.00 0.00 0.028wk 0.00 0.00 0.55 0.00 0.00 0.05 0.07 0.18 0.08 0.42 0.00 0.68 0.83 1.00 1.00 0.98Nor 1.00 1.00 0.37 1.00 1.00 0.92 0.00 0.00 0.00 0.00 0.00 0.15 0.00 0.00 0.00 0.003wk 0.00 0.00 0.05 0.00 0.00 0.05 0.45 0.30 0.73 0.55 0.98 0.10 0.23 0.05 0.02 0.078wk 0.00 0.00 0.58 0.00 0.00 0.03 0.55 0.70 0.27 0.45 0.02 0.75 0.77 0.95 0.98 0.93k2 = 1k2 = 3k2 = 4k2 = 557Figure 2.8: Sammon mapping visualization of histogram profile distances for k1=8. Average pairwise distances were mapped from 30 repetitions of generating histogram profiles from k-medoids clusterings of image histograms. The sixteen injury experiment samples are coloured based on sample condition, and labelled corresponding to Table 2.1. Axes are unit-less dimensions. 58Figure 2.9: Correlation between dispersion and accuracy for CVKNN classification of sample identity. 30 repetitions of CVKNN classification at k1=8, were plotted for k2=1, 3, 4, and 5, along with their linear regression. No significant correlation was found (N=30, R2=.006, .032, .028, and .130 respectively. F-test p-values = .68, .33, .37, and.54 respectively).59Figure 2.10: Accuracy of sample prediction with increasing CVKNN k2 values for selected k1 values. Error bars represent standard error of the mean.60Figure 2.11: Accuracy of CVKNN sample prediction with increasing offset and size penalties. Average accuracy for 60 repetitions of sample prediction at k1=12 and k2=1,3, and 5 using varying Offset Penalty (A), and Size Penalty (B) values. Error bars are stranded error of the mean. Horizontal bars represent significant deviation from the zero penalty condition for each k2 value by Holm-Bonferroni corrected t-test at p=.05 61Figure 2.12: Accuracy of CVKNN sample prediction with varying offset and size penalties for k2=1,3, and 5. k1=12 for each case. Each cell in the heat map is the average of 60 repetitions.62Table 2 Linear model statistics for penalty parameters and their interaction. A linear model was fit to the prediction accuracy when offset penalty and size penalty were varied for k1=12, and k2=1,3, and 5. The estimated coefficient for each penalty, the intercept, and their interaction ( offset penalty × size penalty) along with their error, and student-t statistic comparing the coefficient estimate against an estimate of zero. Coefficient Standard Error tStat pValueIntercept 8.85E-01 5.17E-03 1.71E+02 2.78E-131SizePenalty 4.06E-04 8.75E-05 4.64E+00 1.01E-05OffsetPenalty -3.54E-04 3.23E-04 -1.10E+00 2.76E-01Interaction -1.19E-05 5.46E-06 -2.18E+00 3.12E-02Intercept 8.30E-01 8.14E-03 1.02E+02 1.31E-107SizePenalty 8.70E-04 1.38E-04 6.32E+00 6.28E-09OffsetPenalty -2.63E-03 5.09E-04 -5.17E+00 1.13E-06Interaction 1.42E-05 8.60E-06 1.65E+00 1.02E-01Intercept 7.71E-01 1.17E-02 6.60E+01 6.17E-88SizePenalty 1.49E-03 1.97E-04 7.52E+00 1.82E-11OffsetPenalty -4.63E-03 7.29E-04 -6.34E+00 5.66E-09Interaction 3.13E-05 1.23E-05 2.54E+00 1.25E-02k2=1K2=3K2=52.3.3 Semi-automated thresholding of myelin water spacesFor automating the intensity thresholding of images to separate intact myelin from vacuous water spaces in traced images, images are clustered based on histogram distance using k-medoids and the manual segmentation of the medoid histogram is applied to other histograms in the cluster based on their respective offset from the medoid histogram. Prediction error for each image during model optimization was found by k-fold cross validation.Greedy optimization found the parameters that minimizes total squared prediction error for the injury data-set relative to manual segmentation at k1= 12, Offset penalty = 21, and Size penalty=20. The performance of this model for predicting the training injury data-set is visualized in a correlation and Bland-Altman plot (Figure 2.13). Correlation between this approach and manual thresholding is high (R=0.89, Spearman-Rho p-value= 9.56E-184) and the Bland-Altman plot shows general tendency to underestimate intact myelin fraction (IMF) with a greater magnitude than in cases where it overestimates myelin area which results in an overall methodological bias of 0.0226 IMF and standard error of 0.143 IMF. Overall these data indicate that the model optimization process was able to fit a model to the data-set.Similar performance was found when these parameter settings were used to test the approach with the osmolarity data-set (R=0.88, Spearman-Rho p-value= 2.22E-169, Fisher’s z-transform for correlation comparison with training model p=.491) (Figure 2.14). The Bland-Altman plot showed a similar tendency to underestimate myelin area in a minority of images. Comparing the test model with the training model, the resulting methodological bias was actually lower than in the training case at 0.00836 IMF with a standard error of 0.099 IMF(Table 3). This likely stems from the more homogenous sample conditions of the osmolarity data-set compared to the injury data-set. Overall, this indicates that the optimized parameters are generalizable to similar, but unrelated data-sets. 6364Figure 2.13: Correlation and Bland-Altman plot of model training prediction vs. manual segmentation of intact myelin fraction for semi-automated thresholding. Intact myelin area for each image in the injury data-set was predicted using k-fold cross validation. The resulting predicted intact myelin fraction is plotted against manual segmentation (left) and used to construct a Bland-Altman plot (right). In the correlation the solid line is a linear fit of the data-points, the dashed line indicates an ideal model. In the Bland-Altman plot the average bias of the measurement is reported along with the p-value of a two-tailed t-test for the bias relative to zero, and the standard error of the difference is indicated by dashed lines.65Figure 2.14: Correlation and Bland-Altman plot of model test prediction vs. manual segmentation of intact myelin fraction for semi-automated thresholding. Myelin area for each image in the osmolarity data-set was predicted against the injury data-set model. The resulting predicted intact myelin fraction is plotted against manual segmentation (left) and used to construct a Bland-Altman plot (right). In the correlation the solid line is a linear fit of the data-points, the dashed line indicates an ideal model. In the Bland-Altman plot the average bias of the measurement is reported along with the p-value of a two-tailed t-test for the bias relative to zero, and the standard error of the difference is indicated by dashed lines.Table 3 Comparison of training and test correlation for prediction of semi-automated thresholding.Bias and standard error of prediction for the training and test application of the optimized model are shown along with the p-value for Spearman-Rho rank correlation test, Coefficient of determination, and number of images in the data-setBias Standard error p value nInjury data-set (training) 0.0226 0.143 <.001 0.80 523Osmolarity data-set (test) 0.00836 0.099 <.001 0.78 505R22.3.4 Predicting intact myelin area in unsegmented images based on image histogramsFor automating the prediction of intact myelin fraction based on the intact myelin area of images with similar intensity histograms, images are assigned the average myelin area of the k-nearest images in themodel based on intensity histogram distances. Prediction error for each image during model optimization was found by k-fold cross validation.Greedy optimization of the parameters found the model that minimizes total squared prediction error for the injury data-set relative to manual segmentation at k2= 15, and Offset penalty = 0. The performance of this model for predicting the training injury data-set is visualized in a correlation and Bland-Altman plot (Figure 2.15). Correlation between this approach and manual thresholding is high (R=0.78, Spearman-Rho p-value= 4.84E-10), however the Bland-Altman plot shows a substantial systemic bias: Low myelin area images (myelin area < ~0.2) tend to be overestimated and high myelin area images (myelin area > ~0.5) tend to be underestimated, this is characteristic of a median sample bias for the model. The source of this bias becomes clear when looking at the distribution of manual myelin area in the injury data-set (Figure 2.16). The majority of images are between 0.10 to 0.50 manual IMF and given that the optimized k2 is fifteen, images at the edges of the distribution will necessarily have many 15-nearest neighbours representing median myelin area values if untraced histogram distances does not tightly cluster based on myelin area. Looking at a MDS visualization of untraced injury data-set histogram distances coloured by manual myelin area we can see that there is at best a vague trend towards segregating based on myelin area (Figure 2.17). Having a more even distribution of manual myelin area images in the model may help reduce the systemic bias, however theextent of the overlap means that the standard error for prediction of 0.199 IMF, and correlation coefficient is unlikely to be substantially improved. 66In contrast the overall methodological bias of 0.00265 IMF was not statistically distinct from zero (p=0.548 two-tailed t-test), which is somewhat surprising given the above issues with the model. In spite of said problems the model optimization process was able to fit an accurate, but not precise modelto the data-set which suffers from a median-bias.Performance was substantially worse when the model was tested with the independently sampled osmolarity data-set (R=0.53, Spearman-Rho p-value= 1.16E-37, Fisher’s z-transform for correlation comparison with training model p=1.345E-13) (Figure 2.18). The Bland-Altman plot showed a similar median bias myelin area in prediction. Comparing the test model with the training model, the resulting methodological bias is also higher than in the training case at 0.0455 IMF with a similar standard error of 0.194 IMF(Table 4). This indicates that the optimized model is not generalizable to similar, but unrelated data-sets. In other words, the learned model is overfit to the sample conditions of the injury experiments. 6768Figure 2.15: Correlation and Bland-Altman plot of model training prediction vs. manual segmentation of intact myelin fraction for unsegmented image prediction. Myelin area for each image in the injury data-set was predicted using k-fold cross validation. The resulting predicted intact myelin fraction is plotted against manual segmentation (left) and used to construct a Bland-Altman plot (right). In the correlation the solid line is a linear fit of the data-points, the dashed line indicates an ideal model. In the Bland-Altman plot the average bias of the measurement is reported along with the p-value of a two-tailed t-test for the bias relative to zero, and the standard error of the difference is indicated by dashed lines.Figure 2.16: Histogram of injury data-set images by manually quantified myelin fraction69Figure 2.17: MDS Sammon mapping of untraced injury data-set image histogram distances coloured by manual myelin fraction. Underlines in the figure legend denote inclusive bounds. Axes are unitless dimensions.70Figure 2.18: Correlation and Bland-Altman plot of model test prediction vs. manual segmentation of image myelin fraction for unsegmented image prediction. Myelin area for each image in the osmolarity data-set was predicted against the injury data-set model. The resulting predicted image myelin fraction is plotted against manual segmentation (left) and used to construct a Bland-Altman plot (right). In the correlation the solid line is a linear fit of the data-points, the dashed line indicates an ideal model. In the Bland-Altman plot the average bias of the measurement is reported along with the p-value of a two-tailed t-test for the bias relative to zero, and the standard error of the difference is indicated by dashed lines.Table 4 Comparison of training and test correlation for prediction of unsegmented images by intensity histogram. Bias and standard error of prediction for the training and test application of the optimized model are shown along with the p-value for Spearman-Rhorank correlation test, Coefficient of determination, and number of images in the data-set.Bias Standard error p value nInjury data-set (training) 0.00265 0.199 <.001 0.61 523Osmolarity data-set (test) 0.0455 0.194 <.001 0.28 505R22.3.5 SegNet image segmentation and myelin area predictionFor automating the segmentation of intact myelin area, a SegNet architecture was trained using randomly rotated and shifted images from the injury data-set, and validated against the original images.The trained network was then tested against the osmolarity data-set.At the end of network training prediction accuracy for the training data reached ~76%, with a final validation accuracy of 85%. Most of the improvement in performance occurred during the first 20 epochs of training, with almost no reduction of the error function after 60 epochs. The performance of this model for predicting the training data-set is visualized in a correlation and Bland-Altman plot (Figure 2.19). Correlation between this approach and manual thresholding is high (R=0.79, Spearman-Rho p-value= 2.01E-113) and the Bland-Altman plot shows general tendency to underestimate IMF with a greater magnitude than in cases where it overestimates IMF which results in an overall methodological bias of 0.0389 IMF and standard error of 0.211 IMF. Overall these data indicate that the model optimization process was able to fit a model to the data-set.Performance was substantially worse when the model was tested with the independently sampled osmolarity data-set (R=0.57, Spearman-Rho p-value= 1.24E-45, Fisher’s z-transform for correlation comparison with training model p=1.79E-11) (Figure 2.18). The Bland-Altman plot showed a similar tendency towards greater error when underestimating IMF, and additionally a trend towards overestimating IMF as manual intact myelin area increased. The resulting methodological bias was slightly lower than in the training case at 0.0199 IMF with a similar standard error of 0.231 IMF (Table4). Overall this indicates that the optimized model is not generalizable to similar, but unrelated data-sets.To get a sense of where the SegNet model is making errors we can compare the model’s prediction for intact myelin space with the manually segmented area. In both the training and test data-set the SegNet 71model tended to underestimate intact myelin area in a subset of images which had lower contrast than the majority of images (Figure 2.21A, Figure 2.22A). There were 132 such images in the injury data-setwhich represent two uninjured samples, one 8-week post-injury sample, and one 3-week post-injury sample. Similarly, in the osmolarity data-set there were 131 such images from four samples with phosphate buffer concentrations of 0.0M, 0.023M, 0.047M, and 0.07M (25.2% of training images, 25.9% of testing images.). That this underestimation occurs in low contrast images across a range of sample conditions suggests that the algorithm had difficulty learning how to segment low contrast images in general. It is obvious that features of an image are harder to distinguish as contrast decreases,however this alone is not a satisfactory explanation as the images are still clearly distinguishable from the human perspective, and given that they constitute ~25% of the training data-set their persistent miss-segmentation would contribute a substantial amount of error during training. It seems then that thenetwork did not have enough flexibility during training to learn to distinguish intact myelin in low-contrast images while maintaining accuracy in the rest of the images. In cases where SegNet predicted myelin area accurately there are several common mistakes (Figure 2.21B, Figure 2.22B): Axonal mitochondria and elements of non-axonal processes are often misidentified as myelin, and both positive and negative errors in identifying the exact borders of intact myelin occur frequently. The first mistake is of course indicative of ineffective model training. However, the latter may be a result of more consistent definition of ‘intact myelin’ being applied compared to manual segmentation. Intensity thresholding was used as a shortcut to separate intact myelin from vacuous spaces in the traced myelin area, however uneven illumination, staining, and other imaging artifacts cause the ideal segmentation threshold to vary across an image, so the manual thresholding process requires balancing over-identifying intact myelin in some areas, with under-identifying it in others. Additionally, there is likely minor variability in manual thresholding simply 72because threshold values near the ideal are not perceptively different. The latter source of inconsistencycan be measured by re-thresholding a set of images and comparing the resulting change in intact myelinarea relative to the original thresholding. This was done with a set of 30 random images from the osmolarity data-set. The average absolute error in segmenting threshold intensity was 0.0218 IMF (±0.0035 SEM) out of a maximum intensity of one, this corresponded to an average absolute error of 0.0104 IMF (±0.0019 SEM). Based on these intrinsic sources of error we should expect to see some error in the area of intact myelin simply as a result of the network more consistently applying its definition of intact myelin. An illustrative example of this is circumstances where areas of intact myelinthat were mistakenly missed during manual tracing are ‘mistakenly’ identified as intact myelin by the trained model (Figure 2.23). During training this area of myelin would produce error for the model, however the model could not learn to classify this as not intact myelin because it looks too much like all the other myelin. Thus, this error in manual tracing for the training data was generalized out of the model. A similar outcome should occur with respect to the error in manual thresholding as well although it is less visually striking.In cases where SegNet overestimated myelin area this is usually a result of incorrectly identifying non-axonal elements as intact myelin. In some cases, this occurred when such areas superficially resemble intact myelin: a brighter round interior surrounded by a darker area (Figure 2.21C). In other cases it is not clear how the model is being ‘tricked’ which again points to ineffective learning (Figure 2.21C,Figure 2.22C).7374Figure 2.19: Correlation and Bland-Altman plot of model training prediction vs. manual segmentation of image myelin fraction for SegNet myelin area segmentation. Myelin area for each image in the injury data-set was predicted using the trained neural network. The resulting predicted intact myelin fraction is plotted against manual segmentation (left) and used to construct a Bland-Altman plot (right). In the correlation the solid line is a linear fit of the data-points, the dashed line indicates an ideal model. In the Bland-Altman plot the average bias of the measurement is reported along with the p-value of a two-tailed t-test for the bias relative to zero, and the standard error of the difference is indicated by dashed lines.75Figure 2.20 Correlation and Bland-Altman plot of model test prediction vs. manual segmentation of image myelin fraction for SegNet myelin area segmentation. Myelin area for each image in the osmolarity data-set was predicted using the trained neural network. The resulting predicted image myelin fraction is plotted against manual segmentation (left) and used to construct a Bland-Altman plot (right). In the correlation the solid line is a linear fit of the data-points, the dashed line indicates an ideal model. In the Bland-Altman plot the average bias of the measurement is reported along with the p-value of a two-tailed t-test for the bias relative to zero, and the standard error of the difference is indicated by dashed lines.Table 5Comparison of training and test correlation for prediction of unsegmented images by SegNet. Bias and standard error of prediction for the training and test application of the optimized model are shown along with the p-value for Spearman-Rho rank correlation test, Coefficient of determination, and number of images in the data-set.Bias Standard error p value nInjury data-set (training) 0.0389 0.211 <.001 0.63 523Osmolarity data-set (test) 0.0199 0.231 <.001 0.33 505R276Figure 2.21: Representative images of SegNet performance prediction of intact myelin in the training data-set. Original image on the right, on the left the SegNet prediction is compared with manually determined intact myelin space to overlay correctly predicted intact myelin (yellow), false negative prediction (teal), and false positive prediction (blue). Selected images represent cases where myelin area is underestimated(A), correctly measured(B), and overestimated(C) by SegNet.77Figure 2.22: Representative images of SegNet performance prediction of intact myelin in the test data-set. Original image on the right, on the left the SegNet prediction is compared with manually determined intact myelin space to overlay correctly predicted intact myelin (yellow), false negative prediction (teal), andfalse positive prediction (blue). Selected images represent cases where myelin area is underestimated(A), correctly measured(B), and overestimated(C) by SegNet.78Figure 2.23: Example image of SegNet prediction revealing errors in manual tracing. Original image on the right, on the left the SegNet prediction is compared with manually determined intact myelin space to overlay correctly predicted intact myelin (yellow), false negative prediction (teal), and false positive prediction (blue). Arrow denotes area of myelin which was correctly identified by the model in spite of being incorrectly omitted from manual segmentation.2.3.6 Dense prediction image segmentation and intact myelin area predictionFor automating the segmentation of intact myelin area, a fully convolutional dense prediction architecture was trained using fixed sampled image patches from images of the injury data-set. During training performance was validated using a subset of the constructed patches. After training performance was validated against the original images. The trained network was then tested against the osmolarity data-set.At the end of network training prediction accuracy for the training data reached ~90%, with a final validation accuracy of 88% for the validation patches, indicating a slight trend towards overfitting to the training data. Most of the improvement in performance occurred during the first 12 epochs of training, with almost no reduction of the error function after 20 epochs. The performance of this model for predicting the training data-set is visualized in a correlation and Bland-Altman plot (Figure 2.24). Correlation between this approach and manual thresholding is high (R=0.96, Spearman-Rho p-value= 1.00E-282) and the Bland-Altman plot shows no systemic bias. There is a consistent methodological bias of 0.0493 IMF and standard error of 0.0960 IMF. Overall these data indicate that the model optimization process was able to fit a model to the data-set.Performance was noticeably worse, but still very accurate, when the model was tested with the independently sampled osmolarity data-set (R=0.884, Spearman-Rho p-value= 1.62E-168, Fisher’s z-transform for correlation comparison with training model p=2.22E-16) (Figure 2.18). The Bland-Altman plot showed no systemic bias. Methodological bias was slightly higher than in the training caseat 0.0616 IMF with a similar standard error of 0.0947 IMF (Table 4). Overall this indicates that the optimized model is generalizable to similar, but unrelated data-sets.79Looking at the performance of the dense prediction network compared to manual segmentation we can see that many of the errors which occurred in the SegNet model are substantially improved. Performance in the set of low-contrast images is substantially improved in both the training and test data-set(Figure 2.26A, Figure 2.27A). The area of intact myelin is now broadly correctly defined. Most of the error is in finding the precise border of myelin. In the majority of images axonal mitochondria is no longer incorrectly detected as intact myelin (Figure 2.26B, Figure 2.27B). Finally, dark staining non-axonal elements are less incorrectly identified as intact myelin (Figure 2.26C, Figure 2.27C).8081Figure 2.24: Correlation and Bland-Altman plot of model training prediction vs. manual segmentation of intact myelin fraction for dense prediction myelin area segmentation. Myelin area for each image in the injury data-set was predicted using the trained neural network. The resulting predicted intact myelin fraction is plotted against manual segmentation (left) and used to construct a Bland-Altman plot (right). In the correlation the solid line is a linear fit of the data-points, the dashed line indicates an idealmodel. In the Bland-Altman plot the average bias of the measurement is reported along with the p-value of a two-tailed t-test for the bias relative to zero, and the standard error of the difference is indicated by dashed lines.82Figure 2.25: Correlation and Bland-Altman plot of model test prediction vs. manual segmentation of image myelin fraction for dense prediction myelin area segmentation. Myelin area for each image in the osmolarity data-set was predicted using the trained neural network. The resulting predicted intact myelin fraction is plotted against manual segmentation (left) and used to construct a Bland-Altman plot (right). In the correlation the solid line is a linear fit of the data-points, the dashed line indicates an idealmodel. In the Bland-Altman plot the average bias of the measurement is reported along with the p-value of a two-tailed t-test for the bias relative to zero, and the standard error of the difference is indicated by dashed lines.Table 6: Comparison of training and test correlation for prediction of unsegmented images by dense prediction. Bias and standard error of prediction for the training and test application of the optimized model are shown along with the p-value for Spearman-Rho rank correlation test, Coefficient ofdetermination, and number of images in the data-set.Bias Standard error p value nInjury data-set (training) 0.0493 0.0960 <.001 0.92 523Osmolarity data-set (test) 0.0616 0.0947 <.001 0.78 505R283Figure 2.26: Representative images of dense prediction performance prediction of intact myelin in the training data-set. Original image on the right, on the left the SegNet prediction is compared with manually determined intact myelin space to overlay correctly predicted intact myelin (yellow), false negative prediction (teal), and false positive prediction (blue). Selected images are the same as those in Figure 2.2184Figure 2.27: Representative images of dense prediction performance prediction of intact myelin in the test data-set. Original image on the right, on the left the SegNet prediction is compared with manually determined intact myelin space to overlay correctly predicted intact myelin (yellow), false negative prediction (teal), and false positive prediction (blue). Selected images are the same as those in Figure 2.22.2.3.7 Expected error of machine learning approaches to quantifying sample myelin areaGiven the average standard deviation of manual IMF in all samples from both the injury and osmolaritydata-sets of 0.1043 IMF with a standard deviation of 0.0277 IMF, and the methodological deviation andbias for each machine learning approach from the previous sections it is possible to calculate the expected error in sample average myelin area for each method for given a number of images per sample (Eq.2.3). If we take each method’s bias in the test data-set as representative of how it would bias an unrelated sample while correcting using the training bias then Methodological Bias is the test bias minus training bias. Here this is plotted for each machine learning approach and manual myelin area measurement for samples of average standard deviation, and samples with average standard deviation plus 1.96 standard deviations of the standard deviation (Figure 2.28). Of course, this assumes that methodological deviation is consistent across myelin area values within each method which is not the case for predicting myelin area from whole image intensity histograms due to the median bias, and for the SegNet segmentation of intact myelin due to its inaccuracy in low-contrast images.Expected standard error=Sample SD+Method SD√Number of images+Methodological Bias (21)8586Figure 2.28: Expected myelin area sample error at increasing sample size for each machine learning approach Lines denote expected error for samples of average standard deviation, dots represent the expected error for samples of average plus 1.96 SD standard deviations. 2.3.8 Contrasting features of SegNet and dense prediction architectures.The case could be made that the dense-prediction architecture performed better than SegNet simply because it had more trainable convolutional elements. To investigate this possibility a modified version of the original SegNet architecture was trained on the injury data-set with an increased number of convolutional filters available to the innermost layers of the network (Table 7 2nd column). This resulted in improved training performance and correlation with manual myelin area in the training process and validation to nearly the same level of accuracy. However, in spite of having more trainable convolution weights than the dense prediction network it only showed a moderate improvement in correlation with manual segmentation for predicting intact myelin in the training and testing data-sets. This indicates that increasing the number of high-level filters helped the SegNet network generalize theoverall appearance of myelin. Since the first and last four convolutional layers were unchanged it is unlikely that it accomplished this by using the additional filters to preserve low-level data (Each artificial neuron at this depth is responsible for a 16 by 16 pixels area in the original image.). On inspection the issue with identifying intact myelin in low-contrast images still persisted which suggests that better processing of low-level data are important to make the network robust to differences in sample staining or illumination.One possibility is that the availability of raw low-level texture detail as a result of skip connection is responsible for the dense-prediction network’s improved performance. It is technically possible for SegNet to preserve this information, and this may free the deeper layers of SegNet to focus on higher-level inferences if a skip-connection is added to the network. This was assessed by passing the output of 1st convolutional SegNet layer to a two-output convolutional layer which was fed into the softmax layer along with the regular output of the network (Table 7 3rd column). This allows the network to combine separate predictions by the SegNet, and skip connection halves of the network. The 87regularization of this skip layer was increased by a factor of 10 as was also done for skip connections inthe dense prediction network. During training the modifications as described resulted in an unstable network (loss function oscillated by an order of magnitude between mini-batches), to address this the skip layer had its learning rate reduced by a factor of 100 which allowed training to proceed. The skip connection slightly improved training performance. Correlation with manual myelin area decreased for the training data-set, but increased for the test data-set indicating that the skip-connection helped the network generalize to a minor extent.Another possibility relates to the presentation of data during network training: in the dense prediction network image patches were rotated by multiples of 90°, whereas in SegNet images were rotated by any integer from 0° to 359° (this was done due to hardware limitations). As pixels are arranged in a grid, rotating them by an angle that is not a multiple of 90 necessarily requires interpolation of pixel values which may alter fine details of images. To test whether this impacts network performance SegNet was trained again while restricting training image rotation to multiples of 90° (Table 7 4th column). This improved training, but not validation accuracy, and reduced correlation with manual myelin area quantification in both the training and test data-set indicating that the network overfit to thetraining data. This suggests that any loss of detail from interpolation is either useful for producing a generalized model, or the benefit of augmenting the test data with less constrained rotation more than offsets the cost of interpolation.Based on these experiments it seems that the improved performance of the dense prediction network is related to how its architecture focuses on processing low level details of an image, rather than attempting to build a high level representation as SegNet does. This may in part be due to the nature of the data-sets: there is less structure to TEM images when the tissue in question is degenerating. Additionally, SegNet may simply be unsuitable for image segmentation tasks that rely on fine detail 88due to its imprecise preservation of low level details. The contrast between the performance of the two networks is even more pronounced when we consider differences in the training protocol: SegNet was trained on whole images augmented by rotation, whereas due to hardware constraints the dense prediction network was trained on only 399 patches of each image which were augmented by rotations of multiples of 90° This should have given SegNet an advantage in terms of learning the data, and generalizing to avoid overfitting.8990Table 7: Contrasting SegNet and dense prediction with modified architectures and training protocols. Data from the initial SegNet and Dense prediction architectures are displayed along with several modified SegNet training protocols: a modified SegNet with additional convolutional filters, SegNet with an additional skip connection passing the first layer to the output, and SegNet where the image augmentation during training was restricted to rotate by multiples of 90°. Number of filters and weights are displayed for each layer, with skip connections denoted with an asterisk. Final training and validation accuracy during networktraining is displayed along with the Pearson correlation for the trained network with manual myelin area in the training and test data-set.Initial SegNet Dense PredictionLayer Filters Weights Filters Weights Filters Weights Filters Weights Filters Weights1 64 576 64 576 64* 576 64 576 64* 5762 64 36864 64 36864 64 36864 64 36864 128 737283 64 36864 64 36864 64 36864 64 36864 256* 2949124 64 36864 64 36864 64 36864 64 36864 512 11796485 64 36864 128 73728 64 36864 64 36864 1024 47185926 64 36864 128 147456 64 36864 64 368647 64 36864 256 294912 64 36864 64 368648 64 36864 256 589824 64 36864 64 368649 64 36864 512 1179648 64 36864 64 3686410 64 36864 512 2359296 64 36864 64 3686411 64 36864 512 2359296 64 36864 64 3686412 64 36864 512 2359296 64 36864 64 3686413 64 36864 256 1179648 64 36864 64 3686414 64 36864 256 589824 64 36864 64 3686415 64 36864 128 294912 64 36864 64 3686416 64 36864 128 147456 64 36864 64 3686417 64 36864 64 73728 64 36864 64 3686418 60 36864 64 36864 64 36864 64 3686419 60 36864 64 36864 64 36864 64 3686420 64 36864 64 36864 64 36864 64 36864Total weights 7.01E+05 1.19E+07 7.01E+05 7.01E+05 6.27E+06Training Accuracy ~76% ~84% ~79% ~81% ~90%Validation Accuracy 85% 87% 87% 85% 88%Training Correlation 0.79 0.88 0.77 0.74 0.96Test Correlation 0.57 0.76 0.59 0.47 0.88Modified SegNetSegNet SkipSegNet 90°x RotationChapter 3: ConclusionThe goal of this thesis was to develop and validate the use of a machine learning technique to automate the segmentation of intact myelin robustly in experimentally manipulated tissue samples. Out of the non-parametric approach and the two neural network architectures only the dense prediction neural network produced acceptable results.The general utility of automated image segmentation is to save time, and guarantee consistent image segmentation. Although there is a drop in accuracy compared to manual segmentation this can be compensated for by taking more images. Effectively this approach slightly increases imaging time withthe microscope in exchange for eliminating manual segmentation. The hardware used for training networks in this thesis likely limited the performance of the networks by restricting the number of trainable elements, however this demonstrates that acceptable results are achievable for specific image segmentation tasks without specialized equipment. It is therefore likely that many labs could employ a similar approach to routine image segmentation tasks using their existing data-sets and hardware without the need for investing in additional infrastructure or expertise.3.1.1 Applications of non-parametric image histogram profile techniquesThe non-parametric approach to sample prediction from image histograms may be useful for detecting sample misidentification in image data-sets. Conceptually detecting misidentified samples by transforming the data to a form suitable for clustering is not novel. Similar techniques are used in both bioinformatics to detect mislabelled samples (63), and in computer science to detect miss-classified images (64). The method demonstrated depends on broad changes in intensity levels which would make it unsuitable for experiments where different sample conditions do not produce such a change in images (ie: morphological changes, or changes in density of objects within particular regions of an 91image). It may be possible to develop this technique to detect such changes by first convolving the image with appropriate filters to cause the intensity histogram to capture such changes, but at that pointone may as well simply train a classification network to perform the task.3.1.2 Applications of dense prediction for intact myelin segmentation Neural network machine learning techniques have been applied to other imaging problems in biology, in fact there are several recent papers specifically focusing on segmenting CNS tissue in EM images (58–62). These approaches focus on identifying intact myelin sheaths for the purpose of morphological measurements or cell counting, and were trained on healthy tissue samples. They likely depend on contextual cues specific to healthy intact myelin sheaths and would likely struggle with segmenting myelin in pathological or abnormal states due to disease or experimental manipulation. By contrast the dense prediction approach demonstrated in this thesis is robustly able to measure myelin in non-normalor intact states. Because the data-sets used images taken at a high level of magnification (for the sake ofaccurately capturing the state of degenerating myelin) it could be used to train networks for identifying myelin at any lower level of magnification by shrinking the images to a lower resolution and retraining the network. This could be used as a basis for automating measurements like g-ratio, axon diameter, measuring the area of spared tissue post-injury, or cell classification and counting with further processing. Some of these tasks would require returning to the original micrographs to manually trace additional cellular components of the images (mitochondria, axons, macrophages, oligodendrocyte processes).There are probably still limitations to what can be accomplished by training networks on this data-set given the limited nature of the tissue sampled. To produce a fully generalized myelin TEM image segmentation model, or to go further and automate segmentation of CNS TEM images generally, we would ideally want to combine images from different imaging studies involving a diverse set of 92magnifications, experimental conditions, brain regions, and imaging equipment. This would allow the network to better abstract how myelin and other components of the CNS can appear in different conditions, which would ideally allow it to segment images accurately in any condition. Some of these variables could be addressed computationally. ANNs can be trained to convert images between artistic styles. For example, converting a watercolour painting into an oil painting (69). It may be possible to train a model to convert good quality TEM images into images with low-contrast, weak staining, knife marks, and so on for the purpose of augmenting training data.Such a model could be a useful component for bootstrapping (in the general sense, not the sampling approach) a generalized 3D model for myelin detection for use with emerging volume and 3D electron microscopy technologies. Approaches to producing models for extracting useful morphological information from the massive data-sets produced by these techniques have largely relied upon laborious manual tracing, or by relying on the underlying structure of biological membranes (61,70). If these techniques are to be used in experiments involving pathological conditions the latter approach may not be viable due to the disruption of normal morphology, and the former approach could be avoided by using extant 2D segmentation models to produce training data for 3D models to train on.93Bibliography1. Chen HSM, Holmes N, Liu J, Tetzlaff W, Kozlowski P. Validating myelin water imaging with transmission electron microscopy in a rat spinal cord injury model. NeuroImage. 2017 Jun;153:122–30. 2. Beierlein M. Chapter 17 - Cable Properties and Information Processing in Dendrites. In: Byrne JH, Heidelberger R, Waxham MN, editors. From Molecules to Networks (Third Edition) [Internet]. Boston: Academic Press; 2014 [cited 2023 Feb 10]. p. 509–29. Available from: https://www.sciencedirect.com/science/article/pii/B97801239717910001783. Membrane Potential and the Passive Electrical Properties of the Neuron. In: Principles of Neural Science [Internet]. Fifth Edition. New York, NY: McGraw-Hill Education; 2014 [cited 2023 Feb 19]. Available from: neurology.mhmedical.com/content.aspx?aid=11016775434. Hursh JB. Hursh, J. B. "Conduction velocity and diameter of nerve fibers. Am J Physiol. 1939;127(1):131–9. 5. Sanders FK, Whitteridge D. Conduction velocity and myelin thickness in regenerating nerve fibres.J Physiol. 1946 Sep 18;105(2):152–74. 6. Basic Neurochemistry. 6th ed. Lippincott-Raven; 1999. 7. Aggarwal S, Yurlova L, Simons M. Central nervous system myelin: structure, synthesis and assembly. Trends Cell Biol. 2011 Oct;21(10):585–93. 8. González-Mariscal L, Betanzos A, Nava P, Jaramillo BE. Tight junction proteins. Prog Biophys Mol Biol. 2003 Jan 1;81(1):1–44. 9. Tetzlaff W. The development of a zonula occludens in peripheral myelin of the chick embryo: A freeze-fracture study. Cell Tissue Res [Internet]. 1978 May [cited 2023 Feb 9];189(2). Available from: http://link.springer.com/10.1007/BF0020926910. Peters A. A RADIAL COMPONENT OF CENTRAL MYELIN SHEATHS. J Biophys Biochem Cytol. 1961 Dec 1;11(3):733–5. 11. Gow A, Devaux J. A model of tight junction function in central nervous system myelinated axons. Neuron Glia Biol. 2008 Nov;4(4):307–17. 12. Harauz G, Ishiyama N, Hill CMD, Bates IR, Libich DS, Farès C. Myelin basic protein-diverse conformational states of an intrinsically unstructured protein and its roles in myelin assembly and multiple sclerosis. Micron Oxf Engl 1993. 2004;35(7):503–42. 13. O’Brien JS. Stability of the Myelin Membrane. Science. 1965;147(3662):1099–107. 14. Lazzarini RA, Griffin JW, Lassman H, Nave KA, Miller R, Trapp BD. Myelin Biology and Disorders [Internet]. Elsevier Inc.; 2003 [cited 2017 Oct 6]. Available from: https://cwru.pure.elsevier.com/en/publications/myelin-biology-and-disorders-29415. Marcus J, Popko B. Galactolipids are molecular determinants of myelin development and axo-glial organization. Biochim Biophys Acta. 2002 Dec 19;1573(3):406–13. 16. Campagnoni AT, Skoff RP. The Pathobiology of Myelin Mutants Reveal Novel Biological Functions of the MBP and PLP Genes. Brain Pathol. 2001 Jan 1;11(1):74–91. 17. Myelin | P. Morell | Springer [Internet]. [cited 2017 Oct 6]. Available from: http://www.springer.com/gp/book/978147571832418. Hess MW, Kirschning E, Pfaller K, Debbage PL, Hohenberg H, Klima G. 5000-year-old myelin: uniquely intact in molecular configuration and fine structure. Curr Biol. 1998 Jul 16;8(15):R511–3.19. Velumian AA, Samoilova M, Fehlings MG. Visualization of cytoplasmic diffusion within living myelin sheaths of CNS white matter axons using microinjection of the fluorescent dye Lucifer Yellow. NeuroImage. 2011 May 1;56(1):27–34. 20. Barres BA, Jacobson MD, Schmid R, Sendtner M, Raff MC. Does oligodendrocyte survival dependon axons? Curr Biol CB. 1993 Aug 1;3(8):489–97. 21. Ludwin SK. Oligodendrocyte survival in Wallerian degeneration. Acta Neuropathol (Berl). 1990;80(2):184–91. 22. David S, Miller RH, Patel R, Raff MC. Effects of neonatal transection on glial cell development in the rat optic nerve: evidence that the oligodendrocyte-type 2 astrocyte cell lineage depends on axons for its survival. J Neurocytol. 1984 Dec;13(6):961–74. 23. Tetzlaff W. Tight junction contact events and temporary gap junctions in the sciatic nerve fibres of the chicken during Wallerian degeneration and subsequent regeneration. J Neurocytol. 1982 Oct;11(5):839–58. 24. Zhang K, Jiang M, Fang Y. The Drama of Wallerian Degeneration: The Cast, Crew, and Script. Annu Rev Genet. 2021 Nov 23;55:93–113. 25. Kozlowski P, Raj D, Liu J, Lam C, Yung AC, Tetzlaff W. Characterizing White Matter Damage in Rat Spinal Cord with Quantitative MRI and Histology. J Neurotrauma. 2008 Jun 1;25(6):653–76. 26. Jon Kleinberg and Eva Tardos. Algorithm Design ( 1st Edition) By Jon Kleinberg And Eva Tardos 2005 PDF [Internet]. 2005 [cited 2023 Feb 20]. Available from: http://archive.org/details/AlgorithmDesign1stEditionByJonKleinbergAndEvaTardos2005PDF27. Bishop CM. Pattern Recognition and Machine Learning [Internet]. 1st ed. Springer New York, NY; 2006 [cited 2023 Feb 20]. 738 p. (Information Science and Statistics). Available from: https://link.springer.com/book/978038731073228. Curry HB. The method of steepest descent for non-linear minimization problems. Q Appl Math. 1944;2(3):258–61. 29. Schubert E. HACAM: Hierarchical Agglomerative Clustering Around Medoids -- and its 95Limitations. 30. Lloyd S. Least squares quantization in PCM. IEEE Trans Inf Theory. 1982 Mar;28(2):129–37. 31. Aloise D, Deshpande A, Hansen P, Popat P. NP-hardness of Euclidean sum-of-squares clustering. Mach Learn. 2009 May 1;75(2):245–8. 32. Mahajan M, Nimbhorkar P, Varadarajan K. The Planar k-Means Problem is NP-Hard. In: WALCOM: Algorithms and Computation [Internet]. Springer, Berlin, Heidelberg; 2009 [cited 2017Sep 19]. p. 274–85. (Lecture Notes in Computer Science). Available from: https://link.springer.com/chapter/10.1007/978-3-642-00202-1_2433. Arthur D, Vassilvitskii S. How Slow is the K-means Method? In: Proceedings of the Twenty-second Annual Symposium on Computational Geometry [Internet]. New York, NY, USA: ACM; 2006 [cited 2017 Sep 19]. p. 144–53. (SCG ’06). Available from: http://doi.acm.org/10.1145/1137856.113788034. Arthur D, Manthey B, Röglin H. k-Means has Polynomial Smoothed Complexity. ArXiv09041113 Cs [Internet]. 2009 Apr 7 [cited 2017 Sep 19]; Available from: http://arxiv.org/abs/0904.111335. K-Means++: The Advantages of Careful Seeding [Internet]. ResearchGate. [cited 2017 Sep 19]. Available from: https://www.researchgate.net/publication/220778887_K-Means_The_Advantages_of_Careful_Seeding36. Kaufman L, Rousseeuw PJ. Partitioning Around Medoids (Program PAM). In: Finding Groups in Data [Internet]. John Wiley & Sons, Inc.; 1990 [cited 2017 Sep 19]. p. 68–125. Available from: http://onlinelibrary.wiley.com/doi/10.1002/9780470316801.ch2/summary37. Cheng D, Zhang B. Careful seeding for the k-medoids algorithm with incremental k++ cluster construction [Internet]. arXiv; 2022 [cited 2023 Feb 18]. Available from: http://arxiv.org/abs/2207.0240438. Altman NS. An Introduction to Kernel and Nearest-Neighbor Nonparametric Regression. Am Stat. 1992 Aug 1;46(3):175–85. 39. Goodfellow I, Bengio Y, Courville A. Deep Learning [Internet]. 2016 [cited 2023 Feb 20]. 800 p. Available from: https://mitpress.mit.edu/9780262035613/deep-learning/40. Gidon A, Zolnik TA, Fidzinski P, Bolduan F, Papoutsi A, Poirazi P, et al. Dendritic action potentialsand computation in human layer 2/3 cortical neurons. Science. 2020 Jan 3;367(6473):83–7. 41. Poirazi P, Papoutsi A. Illuminating dendritic function with computational models. Nat Rev Neurosci. 2020 Jun;21(6):303–21. 42. Beniaguev D, Segev I, London M. Single cortical neurons as deep artificial neural networks. Neuron. 2021 Sep 1;109(17):2727-2739.e3. 43. Ioffe S, Szegedy C. Batch Normalization: Accelerating Deep Network Training by Reducing 96Internal Covariate Shift. ArXiv150203167 Cs [Internet]. 2015 Mar 2 [cited 2020 Feb 20]; Availablefrom: http://arxiv.org/abs/1502.0316744. Kohler J, Daneshmand H, Lucchi A, Zhou M, Neymeyr K, Hofmann T. Exponential convergence rates for Batch Normalization: The power of length-direction decoupling in non-convex optimization. ArXiv180510694 Cs Stat [Internet]. 2018 Oct 6 [cited 2020 Feb 20]; Available from: http://arxiv.org/abs/1805.1069445. Santurkar S, Tsipras D, Ilyas A, Madry A. How Does Batch Normalization Help Optimization? ArXiv180511604 Cs Stat [Internet]. 2019 Apr 14 [cited 2020 Feb 20]; Available from: http://arxiv.org/abs/1805.1160446. Jain AK, Prabhakar S, Hong L, Pankanti S. Filterbank-based fingerprint matching. IEEE Trans Image Process. 2000 May;9(5):846–59. 47. Park SI, Smith MJT, Mersereau RM. A new directional filter bank for image analysis and classification. In: 1999 IEEE International Conference on Acoustics, Speech, and Signal ProcessingProceedings ICASSP99 (Cat No99CH36258). 1999. p. 1417–20 vol.3. 48. Dunn D, Higgins WE. Optimal Gabor filters for texture segmentation. IEEE Trans Image Process. 1995 Jul;4(7):947–64. 49. Sklansky J. Image Segmentation and Feature Extraction. IEEE Trans Syst Man Cybern. 1978 Apr;8(4):237–47. 50. Badrinarayanan V, Kendall A, Cipolla R. SegNet: A Deep Convolutional Encoder-Decoder Architecture for Image Segmentation. IEEE Trans Pattern Anal Mach Intell. 2017 Dec;39(12):2481–95. 51. Fully Convolutional Networks for Semantic Segmentation - IEEE Journals & Magazine [Internet]. [cited 2020 Feb 25]. Available from: https://ieeexplore.ieee.org/document/747807252. Sammon JW. A Nonlinear Mapping for Data Structure Analysis. IEEE Trans Comput. 1969 May;C–18(5):401–9. 53. Alonso-Ortiz E, Levesque IR, Pike GB. MRI-based myelin water imaging: A technical review. Magn Reson Med. 2015 Jan;73(1):70–81. 54. MacKay AL, Laule C. Magnetic Resonance of Myelin Water: An in vivo Marker for Myelin. Brain Plast Amst Neth. 2016 Dec 21;2(1):71–91. 55. Kozlowski P, Liu J, Yung AC, Tetzlaff W. High-resolution myelin water measurements in rat spinal cord. Magn Reson Med. 2008 Apr;59(4):796–802. 56. Chen HSM. Validation and optimization of myelin water imaging in a preclinical model of spinal cord injury [Internet]. University of British Columbia; 2016 [cited 2023 Feb 1]. Available from: https://open.library.ubc.ca/soa/cIRcle/collections/ubctheses/24/items/1.03209609757. Chen HSM, Holmes N, Tetzlaff W, Kozlowski P. Effect of Osmolarity on Myelin Water Fraction Measurement in Aldehyde Fixed Spinal Cord Tissue. 58. Janjic P, Petrovski K, Dolgoski B, Smiley J, Zdravkovski P, Pavlovski G, et al. Measurement-oriented deep-learning workflow for improved segmentation of myelin and axons in high-resolution images of human cerebral white matter. J Neurosci Methods. 2019 Oct 1;326:108373. 59. Kreshuk A, Walecki R, Koethe U, Gierthmuehlen M, Plachta D, Genoud C, et al. Automated tracing of myelinated axons and detection of the nodes of Ranvier in serial images of peripheral nerves. J Microsc. 2015 Aug;259(2):143–54. 60. Ilett M, Wills J, Rees P, Sharma S, Micklethwaite S, Brown A, et al. Application of automated electron microscopy imaging and machine learning to characterise and quantify nanoparticle dispersion in aqueous media. J Microsc. 2020 Sep;279(3):177–84. 61. Lee HH, Yaros K, Veraart J, Pathan JL, Liang FX, Kim SG, et al. Along-axon diameter variation and axonal orientation dispersion revealed with 3D electron microscopy: implications for quantifying brain white matter microstructure with histology and diffusion MRI. Brain Struct Funct. 2019 May;224(4):1469–88. 62. Donnat C, Levy A, Poitevin F, Zhong ED, Miolane N. Deep generative modeling for volume reconstruction in cryo-electron microscopy. J Struct Biol. 2022 Dec;214(4):107920. 63. Fort A, Panousis NI, Garieri M, Antonarakis SE, Lappalainen T, Dermitzakis ET, et al. MBV: a method to solve sample mislabeling and detect technical bias in large combined genotype and sequencing assay datasets. Bioinforma Oxf Engl. 2017 Jun 1;33(12):1895–7. 64. Yang Y, Whinston A. Identifying Mislabeled Images in Supervised Learning Utilizing Autoencoder.In 2022 [cited 2023 Jan 26]. p. 266–82. Available from: http://arxiv.org/abs/2011.0366765. Chan CCM, Khodarahmi K, Liu J, Sutherland D, Oschipok LW, Steeves JD, et al. Dose-dependent beneficial and detrimental effects of ROCK inhibitor Y27632 on axonal sprouting and functional recovery after rat spinal cord injury. Exp Neurol. 2005 Dec;196(2):352–64. 66. Otsu N. A Threshold Selection Method from Gray-Level Histograms. IEEE Trans Syst Man Cybern. 1979 Jan;9(1):62–6. 67. Tibshirani R, Walther G, Hastie T. Estimating the number of clusters in a data set via the gap statistic. J R Stat Soc Ser B Stat Methodol. 2001;63(2):411–23. 68. Elkan C. Nearest Neighbor Classification. 2011. 69. Jing Y, Yang Y, Feng Z, Ye J, Yu Y, Song M. Neural Style Transfer: A Review. IEEE Trans Vis Comput Graph. 2020 Nov 1;26(11):3365–85. 70. Cheng HC, Cardone A, Jain S, Krokos E, Narayan K, Subramaniam S, et al. Deep-Learning-Assisted Volume Visualization. IEEE Trans Vis Comput Graph. 2019 Feb;25(2):1378–91. 98AppendixA AppchapA.1 Complexity analysisTo a certain extent a well-constructed algorithm can be idealized as a black-box, where a valid input will always be converted into an appropriate output. However practically there are some limitations to this idea. Computers have limited storage and memory, so they can only handle a finite amount of data at a time, and more importantly most algorithms are only useful if they can produce a result within a limited amount of time. There are algorithms which we know are computationally sound, but are also practically useless because they are unable to produce a result within the lifespan of a human on a largedata input. The primary theoretical approach for analyzing and comparing the efficiency of algorithms is through complexity analysis, which describes how the run time of an algorithm changes in relation tothe size of its inputs. For example, we would say that an algorithm takes n2 time if the number of computations tends to expand on the square of the size of input n. The simplest program which would result in this level of complexity is a loop that iterates over the input elements nested within a loop that iterates over the input elements.A.2 Big O, omega, and theta notationFormal mathematical definitions exist for reasoning about complexity analysis. These are useful for proving things about an algorithm's time complexity, or setting bounds for any algorithm that could be written for a class of problems. For example, any comparison-based sorting algorithm will be at minimum n log(n) time complexity. An algorithm f(x), has an upper bound of g(x), referred to as O(g(x)) (pronounced: “Big-Oh” g(x)), if Eq.(A22) is satisfied (26). Where x is the size of inputs, M is some positive constant, and xo is some particular value of x. To simplify, this equation says that f(x) has an upper bound of g(x) if there is somevalue of x after which g(x) is always greater than or equal to f(x), and you can multiply g(x) by a 99constant M if that helps make f(x) smaller than g(x). f (x)=O (g ( x))as x→∞ iff ∃M > 0 s.t.∣ f (x)∣≤M ∣g( x)∣ for all x≥xo (A22)The definition for the lower bound of an algorithm is similar to the above, except that the inequality between g(x) and f(x) is reversed, and the notation is referred to as Ω(g(x)). The formal definition is included Eq.(A23)(26). The variables here are as in Eq.(A22) and the simplified method of understanding this is similar, except that g(x) must eventually always be lower than f(x) above some x. f (x)=Ω (g ( x))as x→∞ iff ∃M > 0 st∣ f ( x)∣≥M ∣g ( x)∣ for all x≥x (A23)These concepts may be better understood graphically. Consider Figure A2. The algorithm represented by f(x) is perhaps dominated by start-up costs at low input sizes, but after a certain point it begins accelerating upward following a smooth curve as the size of inputs increases. Initially f(x) has a greater execution time than g1(x), however from xo onward g1(x) is always greater than f(x) therefore the time complexity of f(x) has an upper bound of g1(x), or f(x)=O(x2). Likewise, f(x) is always greater than g2(x)therefore the time complexity of f(x) has a lower bound of g2(x), or f(x)=Ω(x). A tight bound for an algorithm, referred to as θ(g(x)), can be found by showing that a particular function g(x) satisfies both Eq. (A22) And Eq.(A23), which is to say that it is valid both as an upper bound and lower bound. Note that xo and M may be changed when checking each bound, and typically M will need to be changed. Consider Figure A1 where the algorithm f(x) from Figure A2 is both greater and less than g(x)=x2 depending on the constant that g(x) is multiplied by. Therefore, the time complexity of f(x) has a tight bound of g(x), or f(x)=θ(x2). 100A.2.1 Practical assessment of time complexityThere are several consequences of the above definitions which can be used to quickly reason about the time complexity of an algorithm in situations where a robust proof is unnecessary. Firstly, there is a hierarchy of time complexities. An algorithm with O(n!) is 'slower' than one with O(2n), which is slower than one with O(nc). Secondly the largest component of an algorithm dominates its complexity. For example, an algorithm where f(x) = x2 + 3x will be O(x2) as the x2 component expands much faster 101Figure A1: Graphical example of theta notation. The input size x is plotted against execution time y for an arbitraryalgorithm f(x). The function y = x2 is plotted at two different scales for comparison.. Scale of axes are arbitrary. Figure A2: Graphical example of Big-O and omega notation. The input size x is plotted against execution time y for an arbitrary algorithm f(x). Two generic functions are plotted for comparison. Scale of axes are arbitrary. than the x component. Finally, as a consequence of the variable M constant multipliers can not change the time complexity of an algorithm. This is to say that O(x2) is equivalent to O(.5x2), so typically any such constants are dropped from the notation.When reading code or pseudo-code for an algorithm the time complexity may be assessed by examining the nature of loops and recursion in the code that varies on the size of inputs. An algorithm with a loop or recursion that iterates over the size of inputs n will be θ(n). Nested loops or recursion have a multiplicative effect on the time complexity, so a loop that iterates over n within a loop that iterates over n will be θ(n2). Operations traversing a binary tree of size n will be θ(log n) in most circumstances (if the tree is completely unbalanced then it would be equivalent to a list making operations on it θ(n)), so filling such a tree would be θ(n log n).The above definitions and examples cover complexity analysis to the extent that it will be used in this work. Typically, this will be brought up when explaining why particular algorithms were used, when comparing different algorithms, or discussing the limits of certain approaches.A.3 Machine LearningBroadly machine learning is an area of computer science concerned with processing data based on a learned model of patterns or inferences an algorithm discovers in the data rather than a set of explicit instructions. For example, a standard alarm clock may be programmed to wake you up at 8:00am in themorning, whereas a machine-learning alarm clock may be programmed to wake you up “so that you arenot late” which may involve processing traffic data, weather, bus routes, and roadwork through some model to determine at what time to set the alarm. The goal of applying machine learning to a problem is to produce a generalized model that can accurately predict the desired output even when a particular set of inputs have not been seen before. To accomplish this the general approach to machine learning involves training the internal workings of a 102model on a set of training data, validating the hyper-parameters of the model to avoid overfitting, and testing its accuracy on an interdependently sampled data-set. There are two broad classes of machine learning algorithms: non-parametric approaches where the training data are used as the model itself and new samples are predicted by relating them to existing samples, and parametric approaches which attempt to generalize a set of parameters from the training data which can be used to interpret new samples to make a prediction. Training a non-parametric model is accomplished simply by the process of collecting samples, while training a parametric model involves adjusting or discovering some numeric parameters which control how the input to a machine learning algorithm is processed to make a prediction. For example, in one-dimensional linear regression: Given a set of n inputs x and corresponding outputs y there is some parameter w can be found to make a prediction yi=wxi that minimizes the error in prediction for some loss function f(w) across the set of samples. A typical loss function might be the sum of square errors which produces the error or loss function to be minimized Eq.(A3)(27). f (w)=12∑i=1n(w xi−yi)2 ∂ f (w)∂w=∑i=1n(w xi−y i) xi (A24, A25)Minimizing the loss function is done with an algorithm called gradient descent which takes advantage of the fact that this function is differentiable with respect to w Eq.(A4)(27), and that a minimum of the derivative function will have a slope of zero. Starting with an arbitrary value for w at iteration t the derivative of the loss function is used to measure the loss function slope of the current model, which is then used to adjust the value of w to produce a model with a better fit for iteration t+1 Eq.(A5)(27,28). This process is repeated until the improvement of loss with each iteration diminishes past some arbitrary threshold. This approach to learning weights will minimize the loss for any loss function that is smoothly differentiable, and convex (the second derivative is always positive).103w t+1=wt−∂ f (w)∂wtf (wt) (A26)Training a model appears fairly straightforward at first, however, there is a fundamental trade-off between how well the model learns the trading data, and how well the model generalizes to unseen data. Models that learn too well will tend to overfit to the particulars of the input data, which is only a sample of a broader distribution of possible input data. For example when attempting to draw a border between two sets of points from overlapping distributions; allowing a more complex border improves prediction of training data, however it does not reflect the underlying distribution which makes prediction of new data from the same distribution less accurate compared to a simpler border(Figure A3). Mitigating overfitting is typically accomplished through hyper-parameters which adjust how specific the model can become. For non-parametric algorithms this usually involves changing how many already-known samples are compared to the new sample to be predicted, as considering more samples typically averages out peculiarities in the data. In parametric algorithms a penalty is applied to the loss function based on the magnitude of the weights based on the assumption that models with larger weights are more likely to be overfitted to the training data. Expanding on the linear regression example above, we may add a w2 expression to the loss function, with an adjustable hyper-parameter λ which allows us to control the strength of this regularization Eq.(A6)(27). This in turn changes the derivative of the loss function which will affect the outcome of the gradient descent optimization Eq.(A7)(27).f (w)=12∑i=1n(w xi−yi)2+ λ2w2 ∂ f (w)∂w=∑i=1n(w xi−y i) xi+ λ w (A27, A28)To select appropriate hyper-parameters part of the training data may be excluded from the model training to be used as a validation data-set. The remaining training data are used to train several models using different hyper-parameters which are then compared by checking their performance when 104predicting the validation data-set. The accuracy of a model in predicting the validation data is an estimate of each model’s ability to predict the test data-set. The robustness of model validation can be improved by a process called cross-validation where the training data are divided into k equal sized groups. Each of the groups is sequentially used as the validation data-set with the k-1 remaining groups used to train the model. The average accuracy of validation across all k groups is used as the validation accuracy for a particular set of hyper-parameters. The most robust form of cross-validation occurs when k=n such that one sample is held back from training at a time then that sample is predicted using the resulting model. This approach provides the best estimate of test-error. But it requires training n models of size n-1 to check each hyper-parameter m which modifies the time-complexity of a machine learning algorithm f(x) by O( mn f(x)) which can significantly increase the processing time for large data-sets or a large number of hyper-parameters. Once the ideal set of hyper-parameters has been found the predictive ability of the model can be tested by measuring the accuracy when predicting the test data-set. The final test does not introduce any additional considerations: however, it is critical that the test data-set does not influence model training, or hyper-parameter selection if it is to accurately measure the accuracy of the model.105Figure A3: Demonstration of overfitting. (Top) Points are sampled from two normal distributions offset by by four units. The ideal prediction boundary is shown as a dashed line. A more complex boundary perfectly separates the distributions shown by a dotted line. (Bottom) Upon resampling the ‘perfect’ boundary makes fourteen errors, while the ideal boundary makes only ten.106A.4 k-means++ initializationThe naive approach to selecting seeds for Lloyd's algorithm is to randomly select k points from among the inputs. However, there are better approaches that attempt to select points that are likely to result in better refinement of cluster definitions, and reduced chances of aberrant results. k-means++ selects seeds to be spread out among the input points (35). Pseudo-code follows:Choose one seed at random from inputs nWhile number of seeds is less than kCalculate Dn as square distance from each n to nearest seed Calculate Dtotal as sum of all DnSelect new seed from n, with each n having probability Dn/Dtotal to be selected.This algorithm naively appears to be O(k2nd) time complexity, which may be slower than standard k-means if i<k. However, all the pairwise distances can be computed on initialization as all the seeds are guaranteed to always be points in n. This allows the initialization to be O(n2d), while seed selection canbe expressed as matrix operations in O(k2n). This makes the entire operation O(n2d) as typically k<√n. Experimentally k-means++ initialized k-means also tends to outperform randomly seeded k-means by reducing the number of iterations i required to converge, and leads to better clusterings(35).107A.5 Backpropagation of error in network trainingA concern you may have noticed with training neural networks is actually computing the gradients for weights in multi-layer networks. The partial gradient equation of the error with respect to one of the weights of the first layer of neurons would quickly become computationally intractable if you attempted to apply the chain-rule of differentiation to each weight through the entire network. To address this an automatic differentiation approach is used whereby error is propagated backwards through the network one layer at a time(27,39). To demonstrate how this is computed we will use the following example network to show how the gradient of a weight for the output layer is computed, and how the gradient of a weight of a previous layer can be computed based on the error computation of theerror ahead of it. For this example, the output layer of the network will be called the final layer to avoidconfusion with the output of a particular neuron (Figure A4).Finding the gradient of w4 with respect to the error of prediction Eq.(A29) can be accomplished using the already provided gradient equations for the prediction error relative to final layer softmax output Eq.(1.18), softmax output relative to final layer net activation Eq.(1.15), and w4 relative to the net activation of the final layer Eq.(1.9).108Figure A4: Sample network for illustrating error backpropagation.∂Error total∂w4=∂ Error final 1∂output final 1∂ output final 1∂active final 1∂active final 1∂w4 (A29)Finding the gradient of w1 with respect to the total error of prediction has a similar appearing equation Eq.(A30). However, the first term can not be solved from an existing equation because the neuron hide3 sends it’s output two the two final layer neurons final1 and final2. The gradient of the total error with respect to the output of hide3 is thus the sum of the gradient of each final layer neuron’s error withrespect to hide3’s output Eq.(A31). ∂Error total∂w1=∂Error total∂outputhide 3∂ outputhide 3∂ activehide 3∂active hide 3∂w1 (A30)∂Error total∂outputhide 3=∂Error final 1∂output hide3+∂Error final 2∂outputhide 3 (A31)Expanding the first term of Eq.(A31) to encompass the three steps of computation occurring in the finallayer gives the following Eq.(A32). The first two terms of this equation are identical to the first two terms already computed in Eq.(A29). Since this is a now a quantity that is used multiple times we will give it a unique identifier j, for the neuron final1 this quantity is jfinal1. The third term of Eq.(A32) the gradient of final1’s net activation with respect to the output of hide3 is just the current value of the weight w4. Applying this logic back into our original equation for the gradient of w1 gives the followingEq.(A33). Where the first term has been replaced by the sum of j multiplied by the weight for each downstream neuron, the second term is the gradient equation of the ReLU with respect to net activationEq.(1.13), and the final term is the gradient of the net activation with respect to w1 Eq.(1.12). Importantly the value of first two terms here is now jhide3 which can be used to calculate the jhide1 and jhide2 for finding the gradient of weights in the preceding layer of neurons. Backpropagation can be viewed as progressively finding the j-jalue for each neuron in the network, then applying the j-value of a neuron to find the gradient of its weights.109∂ Error final 1∂outputhide 3=∂Error final 1∂output final 1∂ output final 1∂ active final 1∂active final 1∂outputhide 3 (A32)∂Error total∂w1=(w4 j final 1+w5 j final 2)∂output hide3∂activehide 3∂active hide3∂w1 (A33)110A.6 Specialized network layersIncluded here is a description of several specialized types of neural network layers used to construct convolutional neural networks which are used to construct networks throughout the paper.Batch Normalization layerThe MATLAB implementation of this layer has two steps: It fits the input data from each incoming filter image to a normal distribution based on each filter image’s mean and standard deviation, The normalized filter images are then each transformed by separate offset and scale weights that are learnedduring training before outputting to the next layer. This is typically placed between a layer of neurons and their corresponding ReLU layer to normalize the output of the neurons. It was originally thought that this would improve training speed by reducing covariate shift (changes to the weights of early layers have an amplified effect on deeper layers, which makes training less stable) (43). However, it has been found that this is not a significant effect of batch normalization, instead it may be improving training by smoothing the landscape of the loss function or by decoupling sign from magnitude when adjusting weights (44,45).Dropout layerThis layer randomly changes some of the input data to a zero during training based on a set probability. This provides a form of regularization which prevents the network from overfitting by relying too much on a subset of the network.Two-dimensional Convolution LayerIn its simplest form this layer applies a series of convolutional filters to an image. A filter is typically a 3-by-3 matrix of weights which are adjusted during training. The input is an image which is convolved by moving the centre of the filter across the image and taking the dot product of the filter and the image111intensity for the overlapping segments of the image and filter. The result of this product is recorded as the filter moves across the image to produce the convolved image (Figure A5). Typically, in artificial neural networks a convolutional layer will apply many separate filters to produce many different convolved images. Additionally, in deeper layers the input may be all the convolved images resulting from a previous convolutional layer which requires filter matrices with an additional dimension of weights for the different layers on the input. For example, the first convolutional layer of the network applied four 3-by-3 convolutional filters to create four filtered images. If the next convolutional layer isintended to produce eight filtered images then it will need eight 3-by-3-by-4 filters, where the extra dimension allows each finter to act on all four of the filtered layers from the previous layer.Max pooling layerA may pooling layer divides a convolved image into pools, typically 2-by-2 non-overlapping squares, and only passes on the highest value in each pool (Figure A6). This reduces the size of the image by a 112Figure A5: Sample convolution of an image with a 2D filter. The filter is centred on the red-bordered pixel of the image, and the result of the convolution calculation for that filter position is filled into the corresponding pixel in the convolved image. factor of two in each dimension, and has a denoising effect on the input image. This means that a 3-by-3 filter size convolutional layer acting on the output of max pooling layer will be acting on a larger arearelative to the initial image. As a convolutional layer typically produces many convolved images with different filters the max pooling layer acts on each convolved image independently.Max Unpooling layerThis layer essentially reverses a max pooling layer. Given the output of a pooled image, the dimensionsof the unpooled image, and a set of indices from a max pooling layer it will map the values in the pooled image to the indices in the unpooled image. This creates a sparsely populated unpooled image (Figure A7). Importantly the pooled image used does not need to be taken from the same layer as the indices. This is to say that you can use a max pooling layer, feed the pooled image through an arbitrary set of layers that end with an image of the same size, and then unpool it based on the initial pooling. 113Figure A6: Sample max-pooling of an image This is useful for maintaining accuracy of coincidence during the pooling phase. For example, suppose we have as part of our network two filters applied to an image: one for horizontal lines, and one for vertical lines. Areas where these occur close together are likely corners which may be a good thing to know for segmenting object edges but an image may have many different objects with corners so some higher level processing incorporating the output of other filters is needed to identify what a corner belongs to, however once we have the higher level information we need to reapply the spatial location of the line filters to find the corner, this can be done with a corner filter that acts on the unpooled image(Figure A8).114Figure A7: Sample unpooling layer115Figure A8: Example of an unpooling layer preserving low-level spacial information for recombination with higher level processing. Squares represent pixels of images and filters in the network. Red lines represent transmission of unpooling indices (red boxes). In the raw images the filter response is rescaled from zero to one to show the contrast in filter response. Normalization rescales the filter output from negative one to one, and ReLU sets values less than one to zero. Note that the unpooled image intensity values emerging from high-level processing are unchanged in this example for the sake of clarity.A.7 Image processing and Convolutional neural networksNow that we have this notion of convolution filters and specialized neural network layers the question arises: why are convolutional filters useful for processing images? And how can these layers be assembled so as to segment images based on a manually segmented data-set.The general utility of convolutional filters is that their output is a reflection of how well a given area resembles the filter itself. For example, a filter with a black upper half and white lower half will have the strongest response to the lower edge of a dark horizontal line, and the weakest response to the upperedge of a dark horizontal line (Figure A9). Early approaches to computer vision involved applying different predefined or learned filters to an image to accomplish general tasks (edge detection, object segmentation) and more particular problems (identifying specific objects, facial recognition)(46–49). Training a segmentation model for a specific task using such an approach may involve defining different objects in terms of their response to a variety of filters. For example, if we want to tell pictures of different cats apart, an image that has a greater response to a vertical line filter is probably a tiger, and one with a greater response to a spotted filter is probably a leopard. However, there are some issues with such an approach: We would need to use filters at several different scales to allow flexibility for cats at different distances, a generic spot filter may not be too useful for distinguishing cheetahs versus leopards, and higher level processing of images where objects are composed of different parts requires some additional model applied on top of the image filters. From these first two concerns we can see that the flexibility and specificity of a filter-based approach depends on the selection of filters. The two-dimensional convolutional layer addresses this by composing the filter from trainable weights. In terms of implementation a convolutional network layer can be imagined as a constrained version of a fully connected network layer where connections outside of the window of the filter are set to zero, and the weights within each position of the filter are the same for each artificial 116neuron(Figure A10)(27,39). This greatly reduces the number of weights compared to a fully connected network which may appear to make it less flexible, but that filters are applied to the whole image means that when they learn something useful during training it is applied everywhere. This helps with generalizing rather than memorizing.With one of these convolutional layers you can imagine that we could train a set of filters that would beable to distinguish between cheetah, tiger, and leopard fur fairly well under ideal circumstances. But wecould imagine a situation where if a leopard is partially obscured by vertical grass the model may get confused as vertical stripes usually mean tigers. The solution to this problem within our one convolutional layer would be to simply add more filters so that there is more flexibility to learn what the different cats look like in different situations and from different orientations. The problem with this approach is that we are not at any point building a higher level representation of any of our cats. We can only predict a particular pixel in an image based on the response of filters centred on that pixel.The way higher-level processing is accomplished within convolutional neural networks is to apply convolutions sequentially. A three-by-three convolutional layer applied after a three-by-three convolutional layer covers a five-by-five area in the original image, and importantly the second layer ofconvolution is performed on the outputs of the first layer’s filters. The second level of convolution mayfor example combine edges defined at the first level into more complex corners or curves. Another level of convolution may allow these curves and corners to be combined to form shapes, and so on. At higher levels when large areas of the image are being processed together we may say that the filters are encoding a sense of objects or the scene in general. A critical component which helps focus convolutional networks towards this sort of hierarchical representation is the use of max pooling layers.Placing a two-by-two pooling layer between two convolutional layers has some desirable effects. First, the input to the second convolutional layer will be a quarter the size of the initial convolution, this 117drastically cuts down on the number of computations as our network gets deeper making the endeavourof training such a network more computationally tractable. Second this quartering will eventually, through successive convolution and pooling get to the point where the area of the filters corresponds to the entire original image. Thus a convolutional network can gradually narrow its output to come to a single prediction based on a hierarchical organization of convolutional filters(Figure A11). This allows us to perform classification of cat images, but the same general convolutional structure is used in imagesegmentation architectures.118119Figure A9: Convolution of an image with a horizontal filter. The six-by-six pixel filter (not shown to scale) is convolved with the sample image on the left. The resulting image is on the right. Brighter areas denote a stronger response to the filter.120Figure A10: Convolution as a constrained fully-connected network layer. A three-by-one convolution is applied toan arbitrary network layer. The inputs to three of the convolutional neurons are shown, with hatching to show inputs with equivalent weights. All other inputs for these neurons fixed at zero, shown as dished lines for neuron B121Figure A11: A simple convolutional neural network architecture for image classification. A series of two-dimensional convolution followed by max pooling is performed with an increasing number of filters until there are 160 neuron activations relating to some general features of the image. A fully connected layer is used to process these features which is processed by a softmax layer for the final prediction.A.8 Multi Dimensional ScalingTo display high-dimensional data on a two-dimensional page requires some transformation of the data into a set of two-dimensional points. One way of accomplishing this is to compute the pairwise distances between points in the high-dimensional space, and attempt to find a set of two-dimensional points that produces a similar table of pairwise distances. This is the general idea of multi-dimensional scaling(MDS). In this paper a particular variant of MDS called Sammon mapping is used which is weighted such that inaccuracies in mapping both large and small distances in the original space are treated the same Eq.(A34)(52). Where the strain for a set of higher dimensional data-points of the matrix X and a particular two-dimensional mapping matrix Z is computed based on comparing distances in the mapping dist(zi, zj), and distances in the high dimensional space dist(xi, xj). Dividing thedifference by the distance in the higher space for each pair of data-points is what normalizes the response of the loss function between large and small distances in the original space. This can be optimized using gradient descent to produce the mapping that minimizes strain starting from a random mapping. The slope of the strain for one lower dimensional point zm along one of the lower dimensions k is shown Eq.(A35)(52).Within the paper this is used to produce figures for high-dimensional data.--strain(Z )= 1∑i=1n∑j=i+1ndist ( xi , x j)∑i=1n∑j=i+1n (dist (xi , x j)−dist ( zi , z j))2dist ( xi , x j)∂ strain(Z )∂ zmk=−2∑i=1n∑j=i+1n dist ( xi , x j)−dist ( zi , z j)dist ( xi , x j)dist (zi , z j)(zmk− z jk )(A34, A35)122
Thesis/Dissertation
eng
Zoology
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Tetzlaff, Wolfram
Application of machine learning techniques for the quantification of myelin ultra-structure in transmission electron micrographs of the rat spinal cord
Text