Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Upload Folder
/
usctheses
/
uscthesesreloadpub_Volume1
/
Adaptive event-driven simulation strategies for accurate and high performance retinal simulation
(USC Thesis Other)
Adaptive event-driven simulation strategies for accurate and high performance retinal simulation
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
ADAPTIVE EVENT-DRIVEN SIMULATION STRATEGIES FOR ACCURATE AND HIGH PERFORMANCE RETINAL SIMULATION by Adi Azar A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (COMPUTER ENGINEERING) May 2012 Copyright 2012 Adi Azar Table of Contents List of Figures vi List of Tables xii List of Algorithms xvi Abstract xvii Chapter 1: Introduction 1 Chapter 2: Related Research 9 2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Simulation Acceleration Strategies . . . . . . . . . . . . . . . . . . . . . 10 2.2.1 Algorithm-Based Acceleration . . . . . . . . . . . . . . . . . . . 11 Distributed Acceleration . . . . . . . . . . . . . . . . . . . . . . 12 Node-Level Acceleration . . . . . . . . . . . . . . . . . . . . . . 12 Event Queue Level Acceleration . . . . . . . . . . . . . . . . . . 13 2.2.2 Hardware-Based Acceleration . . . . . . . . . . . . . . . . . . . 13 2.3 Neural Simulators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 Micro-physiological Simulators . . . . . . . . . . . . . . . . . . 16 CalC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 EONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 PSICS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 VCell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 MCell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.2 Compartmental Simulators . . . . . . . . . . . . . . . . . . . . . 18 GENESIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 NEURON . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 NeuronC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Nodus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 SNNAP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3.3 Network Simulators . . . . . . . . . . . . . . . . . . . . . . . . . 24 SpikeNET . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Brian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 CSIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 EDLUT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3.4 Retinal-Specific Simulations . . . . . . . . . . . . . . . . . . . . 26 Retsim: a Retinal Simulation Using NeuronC . . . . . . . . . . . 26 ii Virtual Retina . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.4 Neural Simulator Comparison . . . . . . . . . . . . . . . . . . . . . . . 28 Chapter 3: Approach 31 3.1 Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2 Selection of a Simulation Strategy . . . . . . . . . . . . . . . . . . . . . 34 3.3 Event-Driven Simulation Acceleration Strategies . . . . . . . . . . . . . 37 3.4 Basic Simulation Acceleration Strategies . . . . . . . . . . . . . . . . . . 40 3.4.1 Lookahead Event-Driven Algorithm . . . . . . . . . . . . . . . . 40 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 The Lookahead Event-Driven Time Bucket . . . . . . . . . . . . 43 3.4.2 Event-Driven Simulation Algorithm With Event Merging . . . . . 44 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 Operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 The Event-Merging Time Bucket . . . . . . . . . . . . . . . . . . 46 3.4.3 Self-Regulating Event-Driven Simulation Algorithm . . . . . . . 46 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4.4 Differences and Similarities Among the Lookahead, Event Merg- ing and Self-Regulating Algorithms . . . . . . . . . . . . . . . . 48 3.5 Established Simulation Acceleration . . . . . . . . . . . . . . . . . . . . 49 3.5.1 Event Data Structure . . . . . . . . . . . . . . . . . . . . . . . . 49 3.5.2 Macro-models For Complex Nodes . . . . . . . . . . . . . . . . 50 3.5.3 Distributed Event-Driven Simulation . . . . . . . . . . . . . . . . 51 3.6 Simulation Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Chapter 4: Baseline Implementation 58 4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.1.1 Inputs and Outputs . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.1.2 Event Creation From the Input . . . . . . . . . . . . . . . . . . . 59 4.1.3 Cellular-Level Simulation . . . . . . . . . . . . . . . . . . . . . 59 4.1.4 Ruby Programming Language . . . . . . . . . . . . . . . . . . . 60 4.2 Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.2.1 Generic Neural Configurations . . . . . . . . . . . . . . . . . . . 61 4.2.2 Retinal Layer Configurations . . . . . . . . . . . . . . . . . . . . 62 4.2.3 Neuron Specific Configurations . . . . . . . . . . . . . . . . . . 63 4.3 Cells to Layer Assignment . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.3.1 Photoreceptor Layer Cell Assignment . . . . . . . . . . . . . . . 64 4.3.2 Subsequent Layer Cell Assignment . . . . . . . . . . . . . . . . . 68 4.4 Baseline Operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.5 Software Architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.5.1 Event-Driven Simulation Strategy . . . . . . . . . . . . . . . . . 71 4.5.2 Neuron Base Class and Extensibility . . . . . . . . . . . . . . . . 72 4.5.3 Event Queue . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 iii 4.6 Built-in Neural Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.6.1 Photoreceptor Cells - Macro-model . . . . . . . . . . . . . . . . 74 4.6.2 Horizontal Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.6.3 Center-Surround Bipolar Cells . . . . . . . . . . . . . . . . . . . 78 4.6.4 Spiking Ganglion Cells . . . . . . . . . . . . . . . . . . . . . . . 78 4.7 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.7.1 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Running the Simulation . . . . . . . . . . . . . . . . . . . . . . . 81 4.7.2 Simulation Output for a Square Image Input . . . . . . . . . . . . 81 4.7.3 Simulation Output For a Letter “D” Image Input . . . . . . . . . . 82 Chapter 5: Experiments 84 5.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.1.1 Retinal Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.1.2 Horizontal Cells Layer Propagation . . . . . . . . . . . . . . . . 84 5.1.3 Simulation Error Computation . . . . . . . . . . . . . . . . . . . 86 5.1.4 Acceleration Algorithm Computation Cost . . . . . . . . . . . . . 86 5.2 Letter “L” Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.2.1 One Frame With Letter “L” . . . . . . . . . . . . . . . . . . . . . 88 Baseline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Lookahead . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Event-Merging . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 Self-Regulating . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.2.2 Ten Frames With Letter “L” With Linear Intensity Increase . . . . 94 Lookahead . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Event-Merging . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 Self-Regulating . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.2.3 Ten Frames With Letter “L” With Non-Linear Intensity Increase . 105 Lookahead . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Event-Merging . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Self-Regulating . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.2.4 Ten Frames Moving Letter “L” With Linear Intensity Increase . . 119 Lookahead . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 Event-Merging . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 Self-Regulating . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 5.2.5 Letter L Performance Trade-off Summary . . . . . . . . . . . . . 138 5.3 Letter “G” Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 5.3.1 Lookahead For Ten Frames With Letter “G” With Linear Inten- sity Increase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 5.3.2 Lookahead For Ten Frames With Letter “G” With Non-Linear Intensity Increase . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5.3.3 Lookahead For Ten Frames Moving Letter “G” With Linear Inten- sity Increase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 5.4 Advanced Inputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 5.4.1 Random Dots . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 iv Self-Regulating . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 Event-Merging . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 5.4.2 Tommy Trojan Head . . . . . . . . . . . . . . . . . . . . . . . . 156 Event-Merging . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 Self-Regulating . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 5.5 Dynamic Time-Bucket Value . . . . . . . . . . . . . . . . . . . . . . . . 164 Chapter 6: Conclusions and Future Work 166 6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 6.2.1 Dynamic Simulation Acceleration Initiatives . . . . . . . . . . . . 168 Dynamic Time Bucket Computation . . . . . . . . . . . . . . . . 168 Intelligent Self-Regulating Algorithm . . . . . . . . . . . . . . . 171 Event Merging Frequency . . . . . . . . . . . . . . . . . . . . . 172 Dynamic Algorithm Selection . . . . . . . . . . . . . . . . . . . 174 6.2.2 Retinal Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 Photoreceptor Cell Micro-Model . . . . . . . . . . . . . . . . . . 175 Bibliography 180 Appendices 192 Appendix A: The Anatomy of the Mammalian Retina: an Overview . . . . . . 193 The Visual Pathway . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193 The Mammalian Retina . . . . . . . . . . . . . . . . . . . . . . . . . . . 195 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195 Retinal Layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198 Retinal Disorders . . . . . . . . . . . . . . . . . . . . . . . . . . 200 Age-Related Macular Degeneration (AMD) . . . . . . . . 201 Retinitis Pigmentosa (RP) . . . . . . . . . . . . . . . . . . 202 Retinal Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202 Photoreceptors . . . . . . . . . . . . . . . . . . . . . . . . . . . 202 Horizontal Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . 204 Bipolar Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205 Amacrine Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . 206 Ganglion Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . 206 Phototransduction in the Dark and Light . . . . . . . . . . . . . . . . . . 207 In the Dark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207 In the Light . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208 Retinal Spatial Encoding . . . . . . . . . . . . . . . . . . . . . . . . . . 208 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209 Appendix B: Retinal Prosthesis Models . . . . . . . . . . . . . . . . . . . . . 211 VLSI Analogs of Neuronal Visual Processing . . . . . . . . . . . . . . . 211 Artificial Silicon Retina (ASR) . . . . . . . . . . . . . . . . . . . . . . . 214 U.S. Department of Energy (DOE) Artificial Retina Project . . . . . . . . 214 v List of Figures 2.1 Moving a time-consuming module to a hardware accelerator . . . . . . . 14 2.2 A neuron represented using a set of compartments [6] . . . . . . . . . . . 19 2.3 One compartment circuit [6] . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 GENESIS object hierarchy [6] . . . . . . . . . . . . . . . . . . . . . . . 20 2.5 NeuronC synapse representation [100]. . . . . . . . . . . . . . . . . . . . 22 3.1 Two connected neurons . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2 Advancing simulation time using time-driven (right), and event-driven (left) simulation strategies. . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3 Baseline event-driven algorithm flowchart . . . . . . . . . . . . . . . . . 39 3.4 Piecewise linear approximation carried out by the lookahead acceleration algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.5 Lookahead event-driven algorithm illustration . . . . . . . . . . . . . . . 42 3.6 Lookahead event-driven algorithm flowchart . . . . . . . . . . . . . . . . 53 3.7 Monotonic response versus non-monotonic (within a time bucket) . . . . 54 3.8 The time bucket size influence on simulation accuracy. . . . . . . . . . . 54 3.9 Event-driven algorithm with event-merging flowchart . . . . . . . . . . . 55 3.10 Event-driven algorithm with event-dropping flowchart . . . . . . . . . . . 56 3.11 Biological retinal photoreceptor response [76] . . . . . . . . . . . . . . . 57 4.1 The cones’ density given the radius from the fovea’s center [55]. . . . . . 65 vi 4.2 The rod’s density given the radius from the fovea’s center [55]. . . . . . . 66 4.3 Cones’ and rods’ distributions [44] . . . . . . . . . . . . . . . . . . . . . 66 4.4 Empty spot matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.5 The fovea’s center cone distribution . . . . . . . . . . . . . . . . . . . . 68 4.6 Cone density of several density rings . . . . . . . . . . . . . . . . . . . . 69 4.7 Photoreceptor layer output given three rings . . . . . . . . . . . . . . . . 69 4.8 Log polar mapping of the photoreceptor cells . . . . . . . . . . . . . . . 70 4.9 Calendar queue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.10 A simple cone macro-model . . . . . . . . . . . . . . . . . . . . . . . . 75 4.11 Cone macro-model adaptation response . . . . . . . . . . . . . . . . . . 76 4.12 Cone macro-model transduction response . . . . . . . . . . . . . . . . . 77 4.13 Ganglion cell center-surround behavior. Image is taken from [84] and modified. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.14 Baseline model (PR: photoreceptor, HC: horizontal cell, BC: bipolar cell, GC: ganglion cell) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.15 Simulation output for a square . . . . . . . . . . . . . . . . . . . . . . . 82 4.16 Horizontal cells output expansion . . . . . . . . . . . . . . . . . . . . . 82 4.17 Simulation output for the capital letter D . . . . . . . . . . . . . . . . . . 83 5.1 Retinal model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.2 Image of the capital letter “L” . . . . . . . . . . . . . . . . . . . . . . . . 89 5.3 Baseline event growth for the one frame with letter “L” input . . . . . . . 90 5.4 Baseline and the lookahead acceleration algorithm queue lengths for the one frame with letter “L” input . . . . . . . . . . . . . . . . . . . . . . . 91 5.5 Baseline and the event-merging acceleration algorithm queue lengths for the one frame with letter “L” input . . . . . . . . . . . . . . . . . . . . . 93 vii 5.6 Baseline and self-regulating algorithm queue lengths for several HC STOP values for the one frame with letter “L” input . . . . . . . . . . . . . . . 94 5.7 Self-regulating algorithm performance for the one frame with letter “L” input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.8 Sequence of letter L images with linear intensity increase (the actual input contains one input image for each frame) . . . . . . . . . . . . . . 95 5.9 Lookahead acceleration algorithm for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.10 The lookahead acceleration algorithm performance for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . 98 5.11 Event-merging queue lengths for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.12 Event-merging performance for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.13 Self-regulating queue lengths at various HC STOP values for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . 103 5.14 Self-regulating algorithm error and speed up for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . . . 105 5.15 Lookahead queue lengths for various photoreceptor time-bucket values for the letter “L” movie with non-linear intensity increase input . . . . . . 108 5.16 The lookahead acceleration algorithm error for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . . . . . . 109 5.17 The lookahead acceleration algorithm performance for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . . . . 110 5.18 The event-merging acceleration algorithm queue lengths for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . 112 5.19 The event-merging acceleration algorithm performance for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . 114 5.20 The self-regulating algorithm performance for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . . . . . . 116 viii 5.21 The self-regulating algorithm performance for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . . . . . . 118 5.22 The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . . . . 120 5.23 The lookahead acceleration algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . . . . 121 5.24 The lookahead acceleration algorithm error for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . . . . . . 122 5.25 The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) . . . . 124 5.26 The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) . . . . 125 5.27 The lookahead acceleration algorithm error for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) . . . . . . 126 5.28 The event-merging acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . 128 5.29 The event-merging acceleration algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . 130 5.30 The event-merging acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) . 132 5.31 The event-merging acceleration algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) . 134 5.32 The self-regulating acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . 135 5.33 The self-regulating algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . . . . . . . . . 137 5.34 Letter G rounded area . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 5.35 The lookahead acceleration algorithm queue lengths for the letter “G” movie with linear intensity increase input . . . . . . . . . . . . . . . . . 143 5.36 The lookahead acceleration algorithm performance for the letter “G” movie with linear intensity increase input . . . . . . . . . . . . . . . . . 144 ix 5.37 The lookahead acceleration algorithm queue lengths for the letter “G” movie with non-linear intensity increase input . . . . . . . . . . . . . . . 146 5.38 The lookahead acceleration algorithm performance for the letter “G” movie with non-linear intensity increase input . . . . . . . . . . . . . . . 147 5.39 The lookahead acceleration algorithm queue lengths for the moving letter “G” movie with linear intensity increase input . . . . . . . . . . . . . . . 149 5.40 The lookahead acceleration algorithm performance for the moving letter “G” movie with linear intensity increase input . . . . . . . . . . . . . . . 150 5.41 The self-regulating acceleration algorithm queue lengths for the random dots input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 5.42 The self-regulating acceleration algorithm performance for the random dots input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 5.43 The event-merging acceleration algorithm queue lengths for the random dots input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 5.44 The event-merging acceleration algorithm performance for the random dots input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 5.45 Tommy Trojan head . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 5.46 Images represent Input, photoreceptor layer output, horizontal-cell layer output and bipolar-cell layer output . . . . . . . . . . . . . . . . . . . . . 157 5.47 The event-merging acceleration algorithm queue lengths for the tommy trojan head input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 5.48 The event-merging acceleration algorithm performance for the tommy trojan head input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 5.49 The self-regulating acceleration algorithm queue lengths for the tommy trojan head input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 5.50 The self-regulating acceleration algorithm performance for the tommy trojan head input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 5.51 A neuron response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 6.1 The rate-of-change factor illustration. . . . . . . . . . . . . . . . . . . . 169 x 6.2 The produced-volume factor illustration . . . . . . . . . . . . . . . . . . 170 6.3 Dynamic time-bucket assignment . . . . . . . . . . . . . . . . . . . . . . 171 6.4 Several connected neurons . . . . . . . . . . . . . . . . . . . . . . . . . 172 6.5 Biological retinal photoreceptor response [76] . . . . . . . . . . . . . . . 175 6.6 Cone connectivity to other retinal cell in the human’s retina (C: cone and H: horizontal cell) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 6.7 Cone processing flowchart . . . . . . . . . . . . . . . . . . . . . . . . . 177 6.8 Michaelis-Menten saturation function [72] . . . . . . . . . . . . . . . . . 178 9 The anatomy of the human eye [36] . . . . . . . . . . . . . . . . . . . . 194 10 The visual pathway [46] . . . . . . . . . . . . . . . . . . . . . . . . . . 196 11 The ten retinal layers [13] . . . . . . . . . . . . . . . . . . . . . . . . . . 198 12 The outer limiting membrane layer position [35] . . . . . . . . . . . . . . 200 13 The figure illustrates an image (left) and the way it is perceived by AMD (center) and RP (right) patients. Images are taken from [88, 63] and merged together. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203 14 Photoreceptors absorption rates. Curve labeled as follows: 498 the absorbance of rods, 420 absorbance of blue cones, 534 absorbance of green cones, and 564 is the absorbance of the red cones [12]. . . . . . . . . . . . . . . 204 15 ON and OFF ganglion cells [50] responding to various center and sur- round stimuli. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209 16 A silicon retina implementation [65] . . . . . . . . . . . . . . . . . . . . 213 xi List of Tables 2.1 Acceleration strategy comparison . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Neural simulator comparison . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.1 ON ganglion cell behavior . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.2 OFF ganglion cell behavior . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.1 Baseline queue lengths with horizontal-cell propagation limitation turned off and on for the one frame with letter “L” input . . . . . . . . . . . . . . 89 5.2 Baseline and lookahead acceleration algorithm queue lengths for the one frame with letter “L” input. QL is the queue length . . . . . . . . . . . . . 91 5.3 Baseline and the event-merging acceleration algorithm queue lengths for the one frame with letter “L” input . . . . . . . . . . . . . . . . . . . . . . 92 5.4 Baseline and self-regulating algorithm queue lengths for several HC STOP values for the one frame with letter “L” input . . . . . . . . . . . . . . . . 92 5.5 Self-regulating algorithm error for the one frame with letter “L” input . . . 93 5.6 Self-regulating algorithm performance for the one frame with letter “L” input 93 5.7 Lookahead queue lengths for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.8 Lookahead algorithm error for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.9 The lookahead acceleration algorithm performance for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . 98 5.10 Event-merging queue lengths for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 xii 5.11 Event-merging performance for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.12 Self-regulating queue lengths at various HC STOP values (SR: self reg- ulating algorithm) for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.13 Self-regulating algorithm error and speed up for the letter “L” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . . . 104 5.14 Lookahead queue lengths for various photoreceptor time-bucket values for the letter “L” movie with non-linear intensity increase input . . . . . . . . 107 5.15 The lookahead acceleration algorithm error for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . 107 5.16 The lookahead acceleration algorithm performance for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . . . . 108 5.17 The event-merging acceleration algorithm queue lengths for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . 111 5.18 The event-merging acceleration algorithm performance for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . 113 5.19 The self-regulating algorithm performance (SR: self-regulating algorithm) for the letter “L” movie with non-linear intensity increase input . . . . . . 115 5.20 The self-regulating algorithm performance for the letter “L” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . 117 5.21 The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . . . . 120 5.22 The lookahead acceleration algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . . . . 121 5.23 The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) . . . . . 123 5.24 The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) . . . . . 123 5.25 The event-merging acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . 127 xiii 5.26 The event-merging acceleration algorithm performance for the moving let- ter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . . . 129 5.27 The event-merging acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) . . 131 5.28 The event-merging acceleration algorithm performance for the moving let- ter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) . . . 133 5.29 The self-regulating acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . 135 5.30 The self-regulating algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) . . . . . . . . . . . 136 5.31 Letter L simulation experiments summary . . . . . . . . . . . . . . . . . . 138 5.32 The lookahead acceleration algorithm performance for a set of experiments 139 5.33 The event-merging acceleration algorithm performance for a set of exper- iments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 5.34 The self-regulating acceleration algorithm performance for a set of exper- iments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 5.35 The lookahead acceleration algorithm queue lengths for the letter “G” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . 142 5.36 The lookahead acceleration algorithm performance for the letter “G” movie with linear intensity increase input . . . . . . . . . . . . . . . . . . . . . . 142 5.37 The lookahead acceleration algorithm queue lengths for the letter “G” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . 145 5.38 The lookahead acceleration algorithm performance for the letter “G” movie with non-linear intensity increase input . . . . . . . . . . . . . . . . . . . 145 5.39 The lookahead acceleration algorithm queue lengths for the moving letter “G” movie with linear intensity increase input . . . . . . . . . . . . . . . . 148 5.40 The lookahead acceleration algorithm performance for the moving letter “G” movie with linear intensity increase input . . . . . . . . . . . . . . . . 148 5.41 The self-regulating acceleration algorithm queue lengths for the random dots input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 xiv 5.42 The self-regulating acceleration algorithm performance for the random dots input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 5.43 The event-merging acceleration algorithm queue lengths for the random dots input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 5.44 The event-merging acceleration algorithm performance for the random dots input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 5.45 The event-merging acceleration algorithm queue lengths for the tommy trojan head input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 5.46 The event-merging acceleration algorithm performance for the tommy tro- jan head input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 5.47 The self-regulating acceleration algorithm queue lengths for the tommy trojan head input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 5.48 The self-regulating acceleration algorithm performance for the tommy tro- jan head input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 6.1 Event merging trade-offs . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 2 The response of both the retinal ON and OFF ganglion cells for various stimuli focused on the center and/or surround of the retinal ganglion cells. . 208 xv List of Algorithms 1 Event-Driven Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2 Time-Driven Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3 Lookahead Event-Driven Simulation . . . . . . . . . . . . . . . . . . . . 42 4 Event-Driven Simulation with Event Merging . . . . . . . . . . . . . . . 45 5 The Baseline Event-Driven Simulation Algorithm . . . . . . . . . . . . . 72 xvi Abstract The biological retina is a complex system with intricate neural processing, dense con- nectivity and massive number of neurons. Unlike most neurons, most retinal neurons respond with graded potentials that will increase the retinal system neural-state simula- tion update rate. The retinal neurons connect densely to each other using various con- nectivity patterns including feedforward, feedback and laterally leading to an even higher neural-state update rate. Considering also the massive number of neurons, modeling the retinal system is very challenging. The lack of medical cures for common retinal disorders has motivated many research groups to study the retina and to attempt to model this complicated system. Some of the retinal projects have created retinal prosthesis devices fabricated on chips, while other groups have created retinal simulators. Retinal simulators are used as powerful learning tools to understand the retinal neurons and their role in vision. For retinal prosthesis research groups, retinal simulators are important tools to determine a suitable quantity of neurons in the prosthetic device design (e.g. number of photoreceptors to allow face recognition). We describe an efficient event-driven retinal simulator that speeds up the simulation by detecting periods of “inconsequential activity” in the neural responses, and carries xvii out approximations or eliminations of simulation events during these periods. We imple- mented several simulation acceleration algorithms including a lookahead event-driven algorithm that carries out the approximations and self-regulating and event-merging event- driven algorithms that eliminate events. We conducted various simulation experiments using a variety of inputs and showed the acceleration and error trade-offs for each input. The performance of each algorithm has been computed relative to the baseline event-driven simulation that does not include any acceleration algorithms. The lookahead event-driven algorithm carries out a piecewise linear approximation to approximate the future response of the node being simulated given its last two responses. In the lookahead event-driven algorithm we propose, events that are scheduled at a node n i with a delivery time less than or equal to t current +TB i can be delivered to n i at t current , such that TB i is n i time bucket. The lookahead carries out an approximation of the time bucket events using the last two computed responses. Reducing computations will increase the simulation speed, but will possibly reduce accuracy. The time bucket value is determined by the modeler. Event merging and self-regulating aim to reduce the number of events in the system to increase the simulation speed while maintaining accurate results. Self-regulating operates to determine if the simulation can drop a low-value event without processing it. A low- value event is an event that does not influence the destination node’s response or the influence is weak. The event-merging acceleration algorithm attempts to create more- complex meta events, each of which entails a set of events that occurs within a time bucket. The algorithm scans the event queue at certain times and attempts to create one event object (meta event) from several events that are scheduled to be delivered to the same destination with delivery times that fall in the same time bucket. xviii The lookahead algorithm performed well for experiments with linearly changing responses. The self-regulating algorithm performed well for slowly changing responses since it reduces the mean number of reported events from the slowly changing nodes. And, the event-merging acceleration algorithm performed well for nodes receiving many events within short periods of time since this increases the probability of merging a larger number of events. Our main contributions are the following: First, we describe a comprehensive outer retinal simulator framework that is flexible and biomimetic, and processes realistic input, producing tangible output. Second, the simulator employs novel event-driven simulation acceleration strategies that are capable of increasing the simulation speed at much higher rate than the error increases. These strategies would be helpful on any continuous system simulation and are not unique to retinal or biological simulations. xix Chapter 1 Introduction The mammalian retina is a light-sensitive neural tissue lining on the back of the eye. It consists of multiple layers of neurons that work together to create multiple representa- tions of the visual scene. The retina’s main tasks are light transduction and visual scene processing. Once light falls onto the eye, the eye focuses it onto the retina. The retinal photoreceptors transduce the light into graded potentials. The transduced action poten- tials will go through three retinal pathways: feedforward, feedback and lateral pathways. The feedforward path is the path of photoreceptors feeding the bipolar cells and bipolar cells feeding the ganglion cells. The feedback path is the path of horizontal cells feeding to the photoreceptors and the path of amacrine cells feeding back to bipolar cells. The lateral pathways are the paths of horizontal cells communicating laterally with all sur- rounding horizontal cells creating a spatially averaged map of potentials proportional to the light intensities in the visual scene, and the path of amacrine cells communicating with other amacrine cells. All the retinal processing including luminance adaptation and edge detection are carried out on those retinal pathways. The lack of medical cures for common retinal disorders has motivated many research groups to study the retina and to attempt to model this complicated system. Some of the retinal projects have attempted to create retinal prosthesis devices fabricated on chips [65, 64, 120, 2, 22, 10, 119, 30], while other groups have created retinal simulators [117, 99]. Fabricated retinal prosthetic devices and retinal simulators both attempt to re-create the functionality of a retinal model to a certain biological extent. The prosthetic chips 1 in general use simplified models of the biological retina due to technology limitations (mainly dissipated power and volume). Some prosthetic projects focus on creating only the photoreceptor output since this is a sufficient cure for Age-Related Macular Degen- eration (AMD) and Retinitis Pigmentosa (RP) disorders [107, 22]. Those two disorders lead to nonfunctional photoreceptors as they lose their ability to transduce light into elec- trical signals. However, there are other prosthetic projects that attempt to model the entire retina electronically but have not been actually implanted in any patients to date [65, 120]. The central nervous system (CNS) is complex, and the retinal neural circuits are no exception. All of the CNS exhibits complex individual neurons, with massive numbers of inputs (synapses) and significant fan-out at the output. Non-linear processing in the dendritic arbor contributes to the complexity. Adding to the complexity of individual neurons is the complexity of the interconnection structure, with feedforward, feedback and lateral connectivity, in addition to the large quantity of synapses, many of which receive unique inputs. In the retina, the lateral connectivity due to horizontal cells and amacrine cells contribute greatly to the retinal complexity. Most retinal neurons output graded potentials, not spikes, adding to the processing complexity. Finally, the scale of the subsystems in the CNS contributes to the complexity. For example, the human retina contains over one hundred million neurons that work simultaneously to process the neural image. Retinal simulators are used as powerful learning tools to understand the retinal neu- rons and their role in vision. For retinal prosthesis research groups, retinal simulators are important tools to determine a suitable quantity of neurons in the prosthetic device design (number of photoreceptors to allow face recognition). Furthermore, simulators 2 can be used to validate a biological retinal model. For visual cortex modelers, the sim- ulator could be a way to generate ganglion cells’ spike trains that are used as inputs to their systems. Most retinal simulators support larger-scale and more-detailed retinal models relative to the prosthetic projects. Some simulators model the various retinal layers mathemat- ically to create each layer’s presumed functionality without modeling any neuron-level simulation [117]. Other simulators use a more-biological approach to mimic the retinal neuronal connections and functionalities at the neuron level [99]. In order to simulate complex retinal structures fully, a retinal simulator must support the complexity of individual neuronal computations, allow the complex interconnection structure required for simulation, and provide adequate performance for retinal simu- lations that encompass enough neurons in the retina to demonstrate the retinal system behavior. Individual neural complexity can be handled by allowing the user flexibility in describing the neural behavior using mathematics and/or modern programming lan- guages. The interconnection structure is difficult to handle, as exemplified by the hori- zontal cell that requires multiple inputs and outputs throughout the span of the cell. In addition, the sheer number of interconnections contributes to the scaling problem, which appears to be the most difficult hurdle in constructing a retinal simulator. Scalability is a fundamental challenge in simulating the biological retina. The biolog- ical retina contains a massive number of complex neurons and the dense interconnection of these neurons complicates the decomposition of the problem into multiple, smaller simulations. In order to simulate an entire retina or significant portions of the retina, more efficient simulations must be possible. It is very important for some simulations 3 to be capable of simulating a significant portion of the biological retina to allow reti- nal modelers to use this platform for realistic experiments. As such, we are willing to sacrifice accuracy to some extent to scale the simulator. Our hypothesis is that it is possible to construct efficient retinal simulations that favor- ably trade simulation scalability for insignificant accuracy loss. The techniques we pro- pose to achieve this modify a conventional event-driven algorithm to reduce the average number of transfer function computations of various neurons and to reduce the aver- age number of events in the simulation by novel strategies for dropping events, merging events and looking ahead for events. These strategies use the sensitivity of each node to its inputs to make simulation acceleration decisions in a way similar to analog circuit simulators that use multiple time-step values [106]. We also use established simulation acceleration strategies that have been used by prominent simulators including calendar queues [109] and macro-models for complex nodes [56]. Finally, we discuss dynamic simulation acceleration to allow the previously described algorithms to adapt efficiently to the various simulation variables. We propose a bio-inspired scalable retinal simulation framework implemented using an event-driven algorithm. We chose an event-driven simulation because, with an effi- cient approach, it can be much faster than a time-driven simulation. However, with a brute-force event-driven simulation for an analog system like the retina, the mean number of events may grow exponentially. Thus, an intelligent event-driven simulation algorithm must be performed to reduce the simulation time significantly from time-driven simula- tion or brute force event-driven simulation. We propose several simulation acceleration strategies and show the trade-offs of simulation speed and accuracy for various inputs. We discuss several event queueing strategies and their memory and search trade-offs. We 4 implemented a “baseline” event-driven simulation of the biological retina using conven- tional techniques, and we computed its performance for various inputs. Then, we applied various simulation acceleration strategies and compared their relative performance to the baseline implementation. We used the Ruby programming language to write the simula- tor code. We call the project Rubtina, which is a play on the words Ruby and retina. The simulation input is a video sequence of images, and the simulation produces the ganglion cells’ optic nerve spike trains and the outputs of various inner retinal layers. Many simulation acceleration strategies have been proposed. Some strategies pro- pose utilizing hardware [108, 14, 15], while others propose algorithmic approaches like parallelism [62, 38, 60]. We propose several novel algorithms that reduce the average number of computations and events in the simulation including lookahead event-driven, event-driven with event merging and self-regulating event-driven. Reducing both the sim- ulation event count and the frequency of transfer function computations will increase the simulation speed but possibly reduce the simulation accuracy. We discuss our approach to measure the performance of each acceleration algorithm relative to the baseline imple- mentation. The lookahead event-driven algorithm carries out a piecewise linear approximation similar to what have been proposed in several analog circuit simulators [115]. In the lookahead event-driven algorithm we propose, events that are scheduled at a node n i with a delivery time less than or equal to t current +TB i can be delivered to n i at t current , such that TB i is n i time bucket. The lookahead carries out an approximation of the time bucket events using the last two computed responses. Reducing computations will increase the simulation speed, but will possibly reduce accuracy. The baseline simulation is a special case of the lookahead algorithm with a time bucket equal to zero. The algorithm accuracy will depend on the length of the time bucket. The bigger the time bucket, the higher the 5 probability of lower accuracy and vice versa. And, the time bucket itself will depend on the node’s response. The algorithm may not be suitable for discrete systems with sparse response, but it is suitable for analog systems with smoothly changing responses, like the retina. We computed the lookahead algorithm performance relative to the baseline for various inputs and discuss the influence of the time bucket’s size on the simulation performance. We describe the time bucket’s dynamic assignment for each node given the node sensitivity to its inputs in the recent simulation history. Furthermore, we describe a queuing mechanism that allows the simulator to efficiently retrieve events that are sched- uled at a certain node within a certain time bucket. We present other event-driven simulation acceleration strategies including event merging, and self-regulating and show their simulation performance trade-offs. These strategies aim to reduce the events in the system to increase the simulation speed while maintaining accurate results. Self-regulating operates to determine if the simulation can drop a low-value event without processing it. A low-value event is an event that does not influence the destination node’s response or the influence is weak. Determining if an event is a low-value event is a key challenge for this algorithm to work efficiently. When a node n i produces an output at time t current that exactly equals the output created at time t current1 , the new event is a low-value event since other nodes already received that event and could buffer that event until n i sends them a new event. However, dropping events can be extended to consider cases in which the output changes a certain percentage as this event may not have an influence on simulation results. The tolerable percentage change of a node’s output that is worth creating an event for depends on the responses of the nodes receiving the event. The algorithm will record the response change for all simu- lation nodes to their inputs using sensitivity tables, and determine the event value given its recent history in the receiving node sensitivity table. The event-merging acceleration 6 algorithm attempts to create more-complex meta events, each of which entails a set of events that occurs within a time bucket. The algorithm scans the event queue at certain times and attempts to create one event object (meta event) from several events that are scheduled to be delivered to the same destination with delivery times that fall in the same time bucket. Rubtina is a system-level simulator in the sense that it helps modelers 1 to understand the functionality of retinal neurons and the interconnected neurons that form layers rather than giving micro-physiological insights. Rubtina is designed to run experiments like the experimentation of the horizontal cells’ roles in luminance adaptation, and demonstra- tion of the sufficient number of photoreceptors to allow face recognition in a prosthetic device. Rubtina is not designed to be a micro-physiological simulator that allows model- ers to determine a cell’s calcium concentrations or understand the neuron’s ion-channel dynamics. Rubtina expects a configuration file that describes the biological retinal model that is provided by the retinal modeler. This model could be based on the human retina, frog retina or a modified version of the human retina, for example. Rubtina is bio-inspired since all the computations are carried out in a way analogous to the biological retina including ganglion-cell firing, center-surround functionality, horizontal-cell spatial aver- aging, and other biological attributes of retinal behavior. Rubtina’s main contributions are the following: Firstly, the simulation employs novel event-driven simulation acceleration strategies that are capable of increasing the simu- lation speed at much higher rate than the error. Secondly, Rubtina is a comprehensive retinal simulator that is flexible, biomimetic, delay and location aware, and processes realistic input, producing tangible output. Rubtina will be described in subsequent chap- ters. 1 A modeler is a Rubtina user who is writing configuration files, and running simulations. 7 The thesis is organized as follows. In Chapter 2, we survey the simulation accel- eration strategies and the neural simulation literature. In Chapter 3, we describe our approach to accelerate Rubtina. In Chapter 4, we describe Rubtina software architec- ture and implementation. In Chapter 5, we demonstrate various simulation experiments for various simulation acceleration algorithms using multiple inputs. In Chapter 6, we present our conclusions and discuss future work. In Appendix 6.2.2, we give a basic overview of the anatomy of the mammalian retina including the retinal system function- ality, retinal disorders, and various retinal cells. In Appendix 6.2.2, we survey prominent retinal prosthesis projects. 8 Chapter 2 Related Research 2.1 Overview Implementing an efficient retinal simulator requires understanding the literature in sev- eral fields including simulation acceleration, neural simulation and understanding the biological retina. In Section 2.2 we survey a few simulation acceleration strategies, in Section 2.3 we survey prominent neural simulators and in Appendix 6.2.2 we describe the biological retina structure and functionality briefly. Simulations are powerful tools to model and analyze systems. Simulations have been used for decades in many fields including traffic engineering [82], reservoir engineering [81], communication networks [110] and electronic circuits [103]. As systems become more complex, modeling and simulating these systems become more complex. Many research groups attempt to model complicated simulations efficiently in various fields using hardware or efficient algorithms or both. In Section 2.2, we discuss various simu- lation acceleration strategies that have been proposed. Neural simulators are classified based on neuron abstraction and scale into three major groups: micro-physiological, compartmental and network simulators [80]. Micro- physiological simulators simulate a single neuron or very few neurons at a high level of detail, while on the other end, network simulators attempt to model a very large num- ber of simple neurons. Micro-physiological simulators are capable of performing very detailed cell dynamics simulations including simulating intracellular calcium dynamics, 9 and a synaptic junction’s molecular mechanisms. Compartmental simulators are capable of simulating individual neurons and small to medium size neural networks using more abstract neuron models relative to the micro-physiological simulators. These simulators are termed “Compartmental” since they use compartment models in simulating the neural system. The compartmental simulators consider each neuron as an electrical cable, and their synapses are electrical connections [45]. The cable is a series of compartments and each compartment is modeled using electrical circuits that consist mainly of resistors and capacitors. Given a network of neurons, compartmental simulators model the network as an electrical circuit and then simulate the circuit itself. In [95], the author compared seven compartmental simulators. Network simulators attempt to model very large num- ber of neurons by simplifying the neuron model and avoiding an iterative neural update for every neuron. There are many simulators in each category. We will discuss some of the most widely used ones in this chapter briefly. Micro-physiological simulators include CalC [68, 16], EONS [101], PSICS [61], VCell [75], and MCell [7]. Compartmental simulators include GENESIS [11], NEURON [17, 43], NeuronC [99], Nodus [26], and SNNAP [34]. Network simulators include SpikeNET [25, 24], EDLUT [91, 54], and Brian [41]. 2.2 Simulation Acceleration Strategies Many research groups have attempted to scale simulations using efficient algorithms or hardware or both. Hardware-based acceleration involves using hardware to increase the simulation speed. Algorithm-based strategies imply algorithmic solutions that aim to reduce the simulation complexity. 10 The simulation design decisions will influence the simulation scalability. Researchers have used various architectures including event-driven and time-driven architectures. Time-driven simulations increment time using a constant time step, then handle all events that occur at the beginning of the time step. The time step has to be selected to be small enough to capture all events, and large enough to allow for scalable simulations. Event- driven simulations increment the simulations time using the next event’s delivery time, which implies ignoring periods of inactivity. Event-driven simulations are suitable for systems with sparse activity, while time-driven simulations are suitable for simulations that involve events with regular inter-arrival times or with continuously changing values (e.g. analog circuits). 2.2.1 Algorithm-Based Acceleration Algorithm-based simulation acceleration strategies can be classified into three major groups: distributed strategies, node-level strategies and event queue optimization. Dis- tributed strategies require more machines, which may increase the simulation cost. Node- level strategies include optimizing the processing node itself to produce output faster and probably with less accurate results. Event queue strategies involve optimizing the event queue to reduce insertion and search complexity. Cadence Virtuoso Spectre is a prominent example of algorithm-based simulation acceleration. Spectre leverages several optimization strategies including efficient use of multi-core machines, multi-rate technology, matrix-solving optimizations and iteration reduction. Virtuoso Spectre claim that these acceleration algorithms provide a consider- able performance boost with accurate SPICE-level simulation [106, 105]. 11 Distributed Acceleration Distributed acceleration strategies involve running the simulation over several machines to increase the simulation speed such that each machine is executing a certain part of the simulation simultaneously. The number of processing nodes that run on each machine varies depending on the simulation design and node complexity. The Blue Brain project used one node (e.g. a node refers to a neuron in a neural simulator) per machine [104], while others in [29] described a distributed system that divides the simulation into several sub-networks, such that each sub-network runs on a separate machine. Distributed simulation must synchronize various parts of the simulation continuously without high I/O cost to achieve speed improvement, and must maintain accurate event temporal order to keep the simulation accurate. The simulation literature is quite rich in various efficient distributed simulation algorithms [62, 73, 20]. Distributed simulations can be carried out on a multi-core machine [96]. Just like any other conventional system, information has to be passed between various cores efficiently to allow various cores to work simultaneously to solve the same problem. Node-Level Acceleration Node-level acceleration strategies involve optimizing the simulation’s nodes to produce results faster. This can be carried out using an alternative macro-model implementation of the node’s transfer function (e.g. the micro-model) [48, 57, 59]. A good macro-model should allow fast computation, and provide accurate results. The node-level acceleration can be done by measuring the node’s sensitivity to the input. Using the sensitivity value, the update frequency of the node can be reduced to a lower rate, which will increase the 12 simulation speed. In Section 6.2.2, we describe a macro-model implementation for the retinal photoreceptors. Event Queue Level Acceleration In large-scale simulations, optimizing the event queue is very important. The event queue must be capable of inserting new events and must fetch the next chronological event effi- ciently. Simulations in the literature have used various data structures including flat arrays [19], priority queues [91, 54], and calendar queues [109]. Any event queue implementa- tion will have trade offs between the complexity of insertion, search complexity and size of memory used. 2.2.2 Hardware-Based Acceleration Hardware-based acceleration implies using hardware to increase the simulation scalabil- ity. Hardware can be used to fully carry out the simulation (emulation), or hardware can be selectively used to carry out a time-consuming part of the simulation as in Figure 2.1. Hardware acceleration is more expensive, but much faster than algorithm-based acceler- ation strategies. Many analog circuit simulation companies have produced commercial simulation acceleration products like Cadence Palladium [15] and Xtreme Server [113]. Cadence Palladium is capable of simulating 256 million gates, which is beyond the capa- bility of any software simulation. Table 2.1 briefly compares hardware and algorithmic acceleration. 13 Figure 2.1: Moving a time-consuming module to a hardware accelerator 14 Aspect/Class Algorithmic Acceleration Partial Hardware Usage Emulation Speed Low High Very High Cost Low High Very High Communication IO None High Low Hardware Mapping Time None High Very High Table 2.1: Acceleration strategy comparison 15 2.3 Neural Simulators 2.3.1 Micro-physiological Simulators Micro-physiological simulators are used to understand the detailed dynamics of the indi- vidual neuron rather than to simulate neural networks. We will briefly describe a few micro-physiological simulators in the next few sections. CalC CalC (Calcium Calculator) is an open-source simulator developed by Victor Matveev using the C++ language [68, 16]. CalC was developed to study the dynamics of calcium diffusion and buffering in the synaptic terminals. CalC simulates the entrance of calcium through channels into a volume, and its diffusion, buffering and binding to the associated receptors. The simulator uses a model defined in a configuration file to create experi- ments. This configuration file is written using a special CalC language that is described in the CalC manual. The configuration files have to be written manually by the modeler since CalC does not offer a GUI. The output of the program can be plotted using Matlab or the Grace 2D plotting tool [42]. The author offers several sample scripts on the CalC website related to synaptic facilitation, calcium diffusion and more. EONS The EONS (Elementary Objects of the Nervous System) modeling platforms were devel- oped at the Laboratory of Neural Engineering of the University of Southern California 16 [101]. The group focus is understanding and modeling the hippocampal synaptic dynam- ics. EONS was designed to be a modeling platform that helps study molecular mecha- nisms that occur at the synaptic junction and the impact of synaptic geometry. EONS has a GUI that consists of a single glutamatergic synapse including both the pre-synaptic and post-synaptic elements. There are two versions of EONS: V1.1, and V2.0. EONS V1.1 suits students who want to study the basic interactions between synaptic elements using a GUI that facilitates learning. In V1.1, all models are predefined so the modeler does not need to write complicated equations to define new models. EONS V2.0 allows the modeler to determine the detailed model and synaptic geometry equations that he desires. PSICS PSICS (Parallel Stochastic Ion Channel Simulator) is an open-source simulator devel- oped by Textensor Ltd using the Java and Fortran programming languages [61]. The project was developed under a contract with Dr Matt Nolan at The Centre for Neuro- science Research at the University of Edinburgh, UK. PSICS is capable of simulating neurons considering the stochastic behavior of ion channels, and considering the detailed geometry of the channels themselves. PSICS reads an input model that describes the sim- ulation behavior and runs the desired calculations. The output of PSICS will be written to a file. The input model consists of several XML files for various input model compo- nents including morphology, and channels. the PSICS team included many experiments on their website including all needed configuration files [112]. VCell VCell (Virtual Cell) was developed by NRCAM (National Resource for Cell Analysis and Modeling) [75] using the Ruby programming language. The modeler creates the 17 simulation using a client from his own computer, and the simulation itself runs over the web on 84 servers with 256 GHz total CPU power and 119 GB total RAM. VCell has been developed to be a tool for a broad range of experiments. The VCell modeler can build complex models with a web-based Ruby interface to specify topology, geometry, molecular characteristics, and relevant interaction parameters. The Virtual Cell automati- cally converts this biological description into a mathematical system of partial differential equations. Then, numerical algorithms are applied to carry out the simulation. MCell MCell (A Monte Carlo Simulator of Cellular Micro-physiology) is a simulator capable of capturing signaling between various cells in the 3D environment [7]. MCell leverages Monte Carlo algorithms to track the molecules as they diffuse through various channels. The MCell modeler must describe his experiment using a special language called MCell Model Description Language (MDL). MCell will carry out the simulation using the pro- vided MDL file, and the experiment output is visualized using a tool like IBM’s Data- Explorer, IRIT, Rayshade, Povray, Pixar’s Renderman, or DReAMM (MCell’s built-in rendering tool). 2.3.2 Compartmental Simulators Compartmental simulators use a compartmental model style to simulate a neural network. They model the neuron as a set of cylindrical pieces (compartments) connected to each other. Figure 2.2 represents a biological neuron and its associated multi-compartment neural model. Each compartment consists of resistors, capacitors, and voltage sources. Figure 2.3 illustrates a compartment circuit. Compartmental simulators offer individual neuron or small to medium size network simulations. 18 Figure 2.2: A neuron represented using a set of compartments [6] GENESIS GENESIS (the GEneral NEural SImulation System) is an open-source simulator written using the C language [11]. It was designed as a general-purpose neural simulator. GENE- SIS is capable of carrying out single cell simulations [8] or small to medium-size network simulations. GENESIS is one of the most popular neural simulators and it is widely used in academia and by many research groups. GENESIS simulations are built from reusable modules that receive inputs from other modules and perform computations on these inputs. These modules can be neurons with a certain morphology or a small neural network. The modeler develops desired neuron models from compartments. These neurons can be connected with each other to form a new module, and this new module can be reused as the modeler desires. GENESIS simulation elements are organized in a tree structure like these in Figure 2.4. The notation 19 Figure 2.3: One compartment circuit [6] is quite similar to the UNIX directories. For example, “/net/cell[2]/dend[1]/GABA” refers to an inhibitory synapse in the first dendrite compartment of the second cell in the network “net”. Figure 2.4: GENESIS object hierarchy [6] There are three ways to define a GENESIS experiment: interactive prompt, script file, and GUI (e.g. XODUS). All the computational modules are the same regardless of the 20 construction method. Using XODUS should be easier for modelers than other simulators since they can visualize what they are building. PGENESIS (Parallel GENESIS) is a parallel form of the GENESIS simulator. It was developed at the Pittsburgh Supercomputing Center to enable the GENESIS neu- ral simulator to run on a variety of parallel machines. It can be used on most Unix or Linux platforms that support PVM (Parallel Virtual Machine) or MPI (Message Passing Interface). PGENESIS could be used for larger simulations to increase the computation speed. NEURON NEURON is an open-source neural simulator written using the C language, it was designed to allow simulations at the neuron level or for small networks [17, 43]. NEU- RON allows modelers to specify the model parameters using neural idioms rather than mathematical formulations, which facilitates simulation experiments that are described using experimental biological data. NEURON modelers can select from several inte- gration methods that are already defined. NEURON experiments can be created using a GUI or using the hoc language (an interpreted language with C-like syntax). NEU- RON researchers have extended hoc libraries to include many various neural modeling functions. NeuronC NeuronC is an open-source neural simulator written using the C language. It offers a computational language that allows simulating neural circuits with biophysical proper- ties including cell morphology, membrane channels and synaptic connections [99]. The NeuronC simulation approach is similar to NEURON in terms of using a variation of the 21 Crank-Nicholson numerical integration method. However, NeuronC allows quite large- scale simulations that include up to 50,000 compartments. The NeuronC modeler writes his script using a computational language, which is similar to the C programming lan- guage. NeuronC parses this script, builds a neural model that represents what is in the configuration file and then carries out the simulation. The script can include a hierarchy of modules to allow module reusability and to reduce the complexity of big simulations. These modules can be built from NeuronC basic elements or from other modules or from both. NeuronC basic elements include sphere, cable, gap junction, synapse and channel. The sphere represents the cell body, the cable represents a dendrite, and the channel is the membrane channel. NeuronC is capable of creating a 3D visual representation of the modeler-defined neural network. The NeuronC compiler’s first task is to model the various modules in the NeuronC script file. This is carried out by translating spheres and cables into compartments. Each sphere will be represented using one compartment while each cable is represented using a series of compartments. A chemical synapse is represented using filters and a transfer function that defines how much neurotransmitter is released as a function of the pre- synaptic voltage. Figure 2.5 represents a NeuronC synapse [100]. Figure 2.5: NeuronC synapse representation [100]. 22 NeuronC differs from the NEURON and GENESIS simulators in a few aspects. Com- pared to NEURON, it is developed for larger scale simulations while NEURON is a better fit for single neuron simulations. NEURON supports the integration of calcium concen- tration for calcium-sensitive channels, which is not a feature offered by NeuronC. Sec- ondly, the numerical method (Gauss-Seidel iteration) used in NeuronC is more capable of simulating more general neural networks that may include loops. Thirdly, the NeuronC library is richer; it contains built-in cells like photoreceptors and offers a 2D stimulus array. Comparing NeuronC to GENESIS, GENESIS runs experiments from configura- tion files while NeuronC uses a programming-language style. Secondly, NeuronC offers a light stimulator while GENESIS does not. The NeuronC distribution includes a script called Retsim that we discuss in Sec- tion 2.3.4. NeuronC libraries include a 2-dimensional light stimulus and a photoreceptor model that facilitate creating retinal and visual experiments including Retsim. The light stimulus includes spots, bars, and sine wave gratings, which are used to excite the pho- toreceptors. The stimulus parameters (location, wavelength, intensity and blur) can be changed. Nodus Nodus is a MAC-based neural simulator developed using the Fortran Language [26]. Nodus is designed to simulate the electrical behavior of a single neuron or small neural networks (up to 200 neurons). The Nodus GUI facilitates creating simulations. SNNAP SNNAP (Simulation of Neural Networks and Action Potentials) is an educational neural simulator developed using the Ruby programming language. It is designed to create 23 simulations using a GUI that does not require any programming experience [34]. The GUI allows modelers to build models, run simulations and visualize the neuron output. SNNAP is capable of simulating realistic models of single neurons using the Hodgkin- Huxley model and small neural networks. The connections among neurons can be made by either electrical or chemical synapses. The chemical synaptic connections are capable of expressing several forms of plasticity including synaptic depression (decrease in amplitude of the postsynaptic potential) and facilitation (increase in amplitude of the postsynaptic potential). 2.3.3 Network Simulators SpikeNET SpikeNET is an open-source spiking neural network simulator developed using the C lan- guage. SpikeNET was designed to carry out very large neural network simulations (mil- lions of neurons) using a limited-biological model of the neuron [25, 24]. The neurons are modeled using the integrate-and-fire model with very limited parameters. SpikeNET neurons communicate instantaneously without any delay. The software uses an event- driven methodology to perform simulations; that is, neurons that are not firing at any time step are inactive. To illustrate how SpikeNET can be used, a multi-scale face recognition network has been described and simulated. The network organization is quite loosely based on the primate visual system. 24 Brian Brian is an open-source spiking neural network simulator written using the python pro- gramming language. The biggest motivation of Brian is to save the modeler time. Both integrate-and-fire models and Hodgkin-Huxley type models can be used in Brian. Con- nectivity between neurons can be randomized and can include delays. Brian simulations are carried out using the Python language. CSIM CSIM is a tool for simulating heterogeneous networks that contain a variety of neural models (leaky-integrate-and-fire neurons, compartmental based neurons, sigmoidal neu- rons) and synapses. It is intended to simulate large-scale networks that may consist of millions of neurons. CSIM is integrated into Matlab, which allows the modeler to create his simulations using the Matlab scripting language. There is not any specific implemen- tation to facilitate retinal simulations. EDLUT EDLUT is an open-source spiking neural network simulator written using the C++ lan- guage [91, 54]. The simulator uses an event-driven approach to reduce the neural update rate. Furthermore, EDLUT creates lookup tables to pre-compute the neural state for the simulation’s neural models to avoid computing the numerical integration online. EDLUT stores the neural events in a priority queue to reduce the complexity of retrieving events. 25 2.3.4 Retinal-Specific Simulations Retsim: a Retinal Simulation Using NeuronC Retsim is a retinal simulation experiment developed using NeuronC. The retinal model is described using a configuration file that lists a table of cells (including photoreceptors, bipolar cells, ganglion cells, horizontal cells and amacrine cells), and their attributes in the table columns and rows respectively. Just like the biological retina, the stimulus is transduced by the photoreceptor. The photoreceptor signals are transmitted to the bipolar cells which in turn transmit their output signals to the ganglion cells (retinal feed forward path). Furthermore, Retsim includes the horizontal and amacrine cells (feed backward paths). Retsim uses an algorithm to create the cells’ morphologies 1 in the retina. The algo- rithm creates these morphologies based on the number of dendrites for each cell, and various branching patterns. Modelers can digitize the cells’ images and create morpholo- gies from them. Programs like NeuroLucida, IJ-MorphDig, and Xfig are used for this purpose. Retsim scripts contain several neural models including a cone-horizontal cell feed- back network, AII amacrine cell network, night vision pathway, and ganglion cell spike generator. The cone-horizontal cell model was developed as an array of 900 cones. There are gap junctions between cones, and chemical synapses between the cones and horizon- tal cells. The various model parameters (synaptic gain and threshold, membrane channel reverse potential) were calibrated manually until light stimulation matched the known physiological data. 1 Cell morphology refers to the form and structure of the cell features 26 Virtual Retina Virtual retina is an open-source retinal simulation program written using the C program- ming language. The program was designed to give visual cortex modelers a way to create a desired retinal ganglion cell output to be used for their cortical experiments. The pro- gram receives a video as an input and produces a series of spikes that represents the ganglion cells’ outputs [117]. The simulator processing is divided into three stages: Outer Plexiform Layer (OPL), Contrast Gain Control, and ganglion layer. The simulator output was compared with the output of the cat X and Y ganglion cells to validate its functionality. These three layers were modeled mathematically using various linear and non-linear filters that are convolved together to create a certain visual functionality. The mathematical parameters can be customized for each simulation for the entire layer level [111]. The simulation is carried out for these filters without considering any neural details. The first stage is a linear filter that simulates the photoreceptors and horizontal cells. Both the center and surround signals are computed and then the output of the layer is computed as a function of both the center and surround signals. The center signal is the photoreceptor output and is defined as a function of light and a series of temporal and spatial filters. The surround signal is the horizontal cell output and is defined as a function of the center signal and a series of spatial and temporal filters. Finally, the layer output is the difference between the center and weighted surround signals. The second stage is responsible for contrast gain control at the level of bipolar cells. This layer receives input from the OPL layer and processes the signals to create contrast sensitivity. The third stage provides additional temporal shaping and creates the ganglion cell output. In this layer, the bipolar signal is rectified and spatio-temporal filtering is applied. The 27 output from rectification and filtering is the input to a set of noisy leaky-integrate-and-fire neurons (nLIF) that create the final spike trains that represent the ganglion cells’ outputs. This layer simulates both phasic and tonic ganglion cells by changing the weight of the transient weight of the temporal filter. 2.4 Neural Simulator Comparison In Table 2.2, we compare relevant simulators for the following aspects: 1 - Synaptic delay support: very important for a retinal simulator as it allow modelers to study the temporal vision in the retina. 2 - Retinal library: facilitates running experiments as modelers can use the library neurons instead of building new neurons. 3 - Phototransduction: the ability of accepting a realistic retinal input (a movie). 4 - Connecting algorithms: offer a systematic way to connect neurons together without requiring the modeler to specify the connections for every neuron. 5 - Model Flexibility: offers modelers the ability to change a neuron at several levels including transfer function, and connections. Some simulators model the retina as a set of layers without considering the internal behavior which will limit any modeler ability to run variety of experiments. 6 - Scalability: it determines the simulator ability to simulate large-scale simulations. Compartmental simulators are usually slower than network simulators as they attempt to simulate neurons at a higher level of details. 28 7 - Modeling approach: refers to the modeling technique used in the simulator. We used the following terminology regarding the modeling approach: NS: network simulator, C: compartmental modeling, M: mathematical formulation of various retinal layers. 8 - Simulator syntax: refers to the scripting language needed to develop simulation experiments. Most simulators expect the modeler to write a script using a pro- gramming language or a syntax specific to the simulator. Some simulators are capable of developing entire experiments using a GUI. 9 - Graded action potential support: refers to the simulator ability to support and scale for analog simulations. Compartmental simulators can simulate analog signals but they cannot scale. While network simulators are not designed to support analog signals. 29 Property/Simulator GENESIS NeuronC SpikeNET Brian Virtual Retina NODUS NEURON EDLUT Rubtina Synaptic Delay Yes Yes No No No Yes Yes Yes Yes Retinal Library No Yes No No No No No No Yes Phototransduction None Limited Image N/A Movie N/A N/A N/A Movie Connecting Algorithms Limited Yes No No N/A N/A N/A No Limited Model Flexibility Yes Yes Yes Yes No Yes Yes Yes Yes Scalability Limited Limited Yes Yes Yes No No Yes Yes Modeling Approach C C NS NS M C C NS NS New Syntax Yes Yes No Python N/A No (GUI) Yes Yes Ruby Graded Potential No No No No No No No No Yes Table 2.2: Neural simulator comparison 30 Chapter 3 Approach 3.1 Statement Given an event-driven simulation for a complex system like the biological retina, the number of scheduled events could grow exponentially over time and scale, challenging the scalability of the simulation, as we will show in the following discussion. Therefore some systematic methods are needed to accelerate the event-driven simulation to produce scalable simulations. The retina’s massive number of neurons, neural complexity and complex connectivity along with graded potentials as outputs are the main reasons that contribute to the astonishing complexity of this system. The human retina contains over one hundred million neurons creating a scalability challenge for retinal modelers. Simulating a large number of nodes can be very important to allow modelers to create meaningful experiments. In the human retina, for example, each ganglion cell 1 receives an input from an average of one hundred photoreceptors. There are nm ganglion cell axons in the optic nerve, to create an output resolution of nm on the ganglion cells, we might need 100nm photoreceptors. Assuming we are going to simulate 1% of the human retina, we needn = m = 100 which requires simulating over 10 6 neurons. Simulating such a large number of neurons to create a resolution of100x100 is 1 The biological retina output is the ganglion cells’ outputs. 31 very challenging, and compartmental and micro-physiological neural simulators are not capable of carrying out such simulations efficiently. The neural complexity in the retina is different from most neural networks since many cells are responding with graded potentials. Unlike most neurons, the majority of retinal cells respond using graded potentials (analog responses), which indicate a high frequency of neural state change 2 that could increase the mean number of scheduled events in the simulation. The input location influence imposes complex retinal models in the simula- tion that will consume more time relative to simple neural models. The retinal cells connect densely to each other using various connectivity patterns including feedforward, feedback and laterally, leading to an increase in the mean num- ber of events in the simulation due to “broadcasting” events in multiple directions and with varying arrival times. Furthermore, the dense neural connections may connect dis- tal neurons together, making it hard to decompose the simulation into several smaller simulations that can run independently. To appreciate the need for event-driven acceleration algorithms, let us consider a reti- nal model R that contains N neurons. Each neuron has an average of S synapses with inputs to the synapses arriving at different times, and produces F events at the input synapses of other neurons with a mean inter-arrival time at theS synapses of seconds as illustrated in Figure 3.1. Neuron i updates its state whenever neuron j changes its state, provided that both neurons are connected. The retinal model R will have NF / updates per second. In the worst-case scenario, the mean number of events grows exponentially with simulation time such that the number of events will increase by a factor ofF every seconds, assuming that an event is produced every time an input to the presynaptic (producing) neuron changes. To illustrate the exponential event growth, let us assume 2 Neural state is the membrane potential. 32 that at timet current every neuron changes its state producingNF events. In the worst- case scenario, every event has a delay different from every other event creating different arrival times for every event. As such, theNF events will arrive at the N neurons at different times leading toNF transfer function computations, which will produceNF 2 events. The newly producedNF 2 events will arrive at the N neurons again at differ- ent times producingNF 3 events, and the event count continues to grow exponentially as the simulation time progresses. Thus, a brute-force event-driven algorithm needs to be enhanced to avoid the exponential growth of events. All retinal neurons except gan- glion cells respond with graded potentials so that potentials in the retina are continuously changing. The rate of input changes drives the entire simulation, so in our simulation we set the rate of sampling inputs, and that determines subsequent event changes in the system. Since the inputs are changing fairly rapidly, without any limitation on input sam- pling rates, that implies ! 0. The low value of leads to frequent time-consuming neural updates in the retinal system R. In practice, we limit the event rate significantly by altering the sensitivity of input changes that trigger events in the simulation. Figure 3.1: Two connected neurons The research approach was determined in large part by choosing an event-driven sim- ulation strategy. Given this choice, a “baseline” simulator has been implemented. Then, 33 several simulation acceleration strategies have been proposed and implemented. The sim- ulation speed and accuracy tradeoffs of the proposed acceleration algorithms have been computed relative to the baseline implementation. Our acceleration strategies approach the problem by reducing the mean number of processed and scheduled simulation events intelligently to favorably increase simulation speed with insignificant accuracy loss using several novel algorithms including lookahead event-driven, event-driven with event merg- ing and self-regulating event-driven. Then, we used several established simulation accel- eration strategies including macro-models, and an efficient event queue. Finally, we pro- posed dynamic simulation acceleration strategies that allow the acceleration algorithms to dynamically adapt to the simulation efficiently. Our contributions are the following: Firstly, the simulation employs novel event- driven simulation acceleration strategies that are capable of reducing the average number of events in the simulation without sacrificing significant accuracy. Secondly Rubtina is a comprehensive retinal simulator that is flexible, biomimetic, delay and location aware, and processes realistic inputs, producing tangible outputs. 3.2 Selection of a Simulation Strategy Simulators can be classified using their approach of advancing simulation time into time- driven and event-driven simulators. Advancing simulation time is directly related to the simulation speed and accuracy since simulators process simulation nodes once they advance simulation time. In time-driven simulators, the simulation time is advanced using a fixed period of time called the time step, and simulators typically process every single node in the simulation 34 every time step. Time-driven simulators choose small time steps to increase the simula- tion accuracy, and the time step value needs to be large enough so that the simulation does not spend too much time carrying out computations for every node in the simulation. The time-driven simulation strategy is usually used for analog system simulations, as inputs and output are continuously changing. The disadvantage of the time-driven simulation strategy is that the simulation will process every node every time step regardless of the node’s activity. In event-driven simulators, the simulation time is advanced using the next event deliv- ery time. This approach allows the simulator to skip periods of inactivity to avoid spend- ing simulation time processing inactive nodes. An event-driven strategy is usually used for discrete systems since these systems have periods of inactivity. Figure 3.2 represents the way time-driven and event-driven simulation strategies advance the simulation time. Algorithms 1 and 2 represent simple algorithms of both the event-driven and the time-driven simulation algorithms respectively. Figure 3.2: Advancing simulation time using time-driven (right), and event-driven (left) simulation strategies. 35 Algorithm 1 Event-Driven Simulation 1: Load input events into the queue 2: initialize simulation clock to zero 3: while event queue is not empty do 4: define current event as the event with the earliest delivery time 5: carry out a transfer function computation to compute the new state 6: generate new events to notify connected nodes about the new state 7: add new events to the event queue 8: set simulation clock to the current event delivery time 9: delete current event from the queue 10: end while Algorithm 2 Time-Driven Simulation 1: initialize the time step 2: initialize simulation clock to zero 3: while simulation clock is less than the simulation entire time do 4: carry out the integral of the transfer function on the period [simulation clock, sim- ulation clock+time step] 5: increment the simulation clock by time step time units 6: end while 36 We have chosen an event-driven simulation strategy for the following reasons: 1. Time-driven simulators don’t scale well because they simulate every node at every time step. 2. Any acceleration that can be applied to a time-driven simulator, can be also applied to an event-driven simulator. 3. Event-driven simulators scale since they skip periods of inactivity, while time- driven simulators do not. Having analog responses in the majority of the retinal cells indicates that there are not frequent periods of inactivity in the simulation. However, there are periods of “inconse- quential activity”, that is, activity that does not make a difference in the simulation that worth even spending simulation time on. In Section 3.3, we propose several algorithms that attempt to identify periods of inconsequential activity in the simulation and compress or remove inconsequential events to create an efficient event-driven simulation algorithm. The inconsequential activity period response computation for a particular node may be fully ignored or approximated using a linear transfer function. 3.3 Event-Driven Simulation Acceleration Strategies We implemented several simulation acceleration strategies for the baseline event-driven simulation that will scale the simulation efficiently as simulation size increases without significant accuracy loss. Figure 3.3 represents the flowchart of the baseline event-driven algorithm that we need to scale. We divide the simulation acceleration strategies into the following: 1. basic simulation acceleration strategies, 37 2. established simulation acceleration strategies, and 3. dynamic simulation acceleration strategies. The basic simulation acceleration strategies attempt to identify periods of inconse- quential activity in the analog responses. These algorithms include the lookahead event- driven, event-driven with event merging and self-regulating event-driven. These three algorithms identify the periods of inconsequential activity, and attempt to reduce the simulation activity for these periods by approximating the node response for the incon- sequential activity periods, merging the inconsequential activity period events together and ignoring low-value events that may not impact the simulation. These algorithms will reduce the mean number of computations carried out by simulation nodes, and reduce the mean queue length. These algorithms will be described in Section 3.4. The established simulation acceleration strategies include simulation acceleration strategies that have been used by prominent simulators. We will use an event queue similar to the calendar queue that has been used by the network simulator NS2 [109]. The simulation literature is rich in distributed and parallel simulation techniques [60, 29], and we discuss the challenges of distributing an event-driven simulation into several processors to carry out a significant retinal simulation. We will use macro-models for retinal complex neurons to reduce the computation time [56]. These algorithms will be described in Section 3.5. Dynamic simulation acceleration strategies will attempt to dynamically modify the simulation behavior to increase the simulation performance. These strategies will be discussed in Section 6.2.1. Neural compartmental simulators like GENESIS [11], NEURON [17, 43], and Neu- ronC [99] compute the neural state for every neuron at every time step. The iterative 38 Figure 3.3: Baseline event-driven algorithm flowchart update for every single neuron reduces the scalability of these simulators. Some of these simulators evolved to include parallel packages like PGENESIS to distribute simulations over several machines [85]. However, our research goal is to scale the simulation itself efficiently and parallel/distributed processing can be added to our efficient simulator. Net- work simulators like SpikeNET [25, 24], and EDLUT [91, 54] propose using an event- driven approach for spiking neural networks. The event-driven approach is motivated by the sparse spiking behavior in the neural system. However, these simulators were not designed to support scalable simulations for systems with graded potential responses like the retina. Many simulators have used an event-driven strategy including the network simula- tor NS2 [110]. However, using a baseline event-driven simulation for a complex system like the retina with graded potentials and complex connectivity will not suffice to scale. Hence, we have added basic simulation acceleration strategies that identify the periods 39 of inconsequential activity and carry out response approximations or event merging or event dropping, which will reduce the simulation event growth. Many time-driven simulators have been proposed using multiple time-steps, includ- ing Cadence Spectre [105, 106], to reduce unnecessary computations. The approach can simply classify each electronic circuit element into one of two groups of elements, fast elements with a small time step and slow elements with a big time step. The basic simula- tion acceleration strategies we propose in Section 3.4 are quite analogous to the multiple time-step approach such that each node’s sensitivity to its inputs is quantified and used to determine possible simulation acceleration, which will reduce the mean number of transfer function computations and the queue mean size respectively. 3.4 Basic Simulation Acceleration Strategies 3.4.1 Lookahead Event-Driven Algorithm Overview Whenever a node receives an event in an event-driven simulation, it carries out a new state computation [115]. As the node’s transfer function complexity increases, computing its new state becomes a time-consuming operation. In large-scale simulations, a node is continuously receiving events and continuously updating its state. The lookahead algo- rithm attempts to reduce the frequency of a node’s computation by approximating the responses of incoming batches of events using the previously computed node responses. Each event batch contains events with delivery times that fall within a time bucket. When a node receives this batch of events, it approximates its responses at the arrival time for every event that falls within the time bucket using the piecewise linear function. Thus, 40 all the events within this time bucket are approximated, and there is no need to carry out a transfer function computation for each individual event. The piecewise linear approx- imation for analog system simulation has been used in analog circuit simulations [115], and the lookahead event-driven algorithm is Rubtina’s algorithm for the linear piecewise approximation for event-driven simulations. Figure 3.4 represents the piecewise linear approximation that is carried out by the lookahead event-driven simulation algorithm. Figure 3.4: Piecewise linear approximation carried out by the lookahead acceleration algorithm Figure 3.5 illustrates the operation of the algorithm. Given a nodeN 7 that is receiving the eventE 37 in the baseline event-driven simulation algorithm. In the lookahead event- driven simulation algorithm, N 7 will get several events falling within its time bucket. Thus, N 7 will receive the set of eventsfE 37 ;E 57 ;E 47 g since all the delivery times for every event in this set fall within t current +TB 7 such that TB 7 is the node N 7 ’s time bucket. Notice that event E 17 was not delivered to N 7 , as its delivery time did not fall withinN 7 ’s current time bucket. The lookahead event-driven simulation starts with a time bucket equals to zero, behaving exactly like baseline event-driven simulation when the time bucket is zero. The 41 Figure 3.5: Lookahead event-driven algorithm illustration simulation’s time is advanced to the batch’s earliest event delivery time. Algorithm 3 and Figure 3.6 represent the operation of the lookahead event-driven algorithm. Algorithm 3 Lookahead Event-Driven Simulation 1: Load input events into the queue 2: initialize simulation clock to zero 3: initialize empty approximations table 4: while event queue is not empty do 5: define current event as the event with the earliest delivery time 6: define current neuron as the current event destination 7: if there is an approximation for the current event then 8: read the event approximation (computed already using the lookahead algorithm) 9: generate new events to notify post-synaptic neurons about the new state 10: else 11: define previous state as the neuron current state 12: carry out a transfer function computation to compute the new state 13: define current state as the neuron’s current state after the recent state update 14: generate new events to notify post-synaptic neurons about the new state 15: define time bucket events as the set of events that are delivered to the current neuron within the next time bucket seconds 16: carry out piecewise-linear approximation. The algorithm uses the previous state and current state to create the linear equation 17: add the approximations of the event responses to the approximation table 18: end if 19: set simulation clock to the current event delivery time 20: delete current event from the queue 21: add new events to the event queue 22: end while Assumptions The lookahead event-driven algorithm will function properly given the following assump- tions: 42 1. The node’s response is monotonic within the time bucket. Non-monotonic response may cause significant accuracy loss. Figure 3.7 illustrates two responses, mono- tonic and non-monotonic. 2. Given a noden i that receives events from a set of nodesfn 0 ;n 1 ;::;n m g, we must satisfy: TB i <minfnd 0 ;nd 1 ;::;nd m g. Such thatnd m is the delay for noden m ’s output to reach noden i , andTB i isn i ’s time bucket. This constraint ensures that events are processed in order. The Lookahead Event-Driven Time Bucket The time bucket is the period of time for which we can execute a batch of events for a certain node in the near simulation future without sacrificing significant simulation accuracy. The time-bucket value will influence the simulation speed and accuracy. With a time bucket equals to zero, the algorithm is the baseline event-driven simulation. A large time bucket will increase the simulation speed but may decrease the accuracy. The simulation speed is increased by a factor ofTB= such that is the event inter- arrival time and TB is the simulation’s average time bucket. The error upper bound can be computed by assuming that all events within a time bucket will produce the output value of the event with the latest delivery time. Then, the error upper bound will be equal to the integral of the produced area (triangle) divided by the time-bucket value. Considering the two points (X 0 ;Y 0 ) and (X 1 ;Y 1 ), the error upper bound will be: 1 2 (Y 1 Y 0 ). The value of(Y 1 Y 0 ) is directly proportional to the time-bucket value. Thus, the lookahead event-driven simulation algorithm speed up and error upper bound are proportional to the time bucket and will increase linearly with the time bucket size increase. Fortunately, the 43 error upper bound is quite pessimistic, so speed up in practice increases faster than error. Figure 3.8 represents two different time buckets used for the same node response. The bigger the time bucket, the faster the simulation and the bigger the error. 3.4.2 Event-Driven Simulation Algorithm With Event Merging Overview The lookahead event-driven simulation algorithm increases the simulation speed by reducing the computation rate. However, the algorithm does not reduce the search and insertion complexity of the event queue. Event merging is designed to reduce the com- plexity of the event queue by reducing the mean number of events using a methodology similar to the lookahead event-driven simulation algorithm. Event merging implies merging a set of events to create a bigger meta event such that this bigger meta event can represent the set of events accurately. The meta event data structure is quite different from the original event as it reports several inputs instead of a single input in the original event. The algorithm merges m events together if they are being sent to the same destination node and their delivery times are within a certain time bucket. The new-formed meta event will be inserted into the queue instead of inserting the original m events, reducing the mean queue length, which will reduce the search and insertion complexity. Operation The event-driven algorithm with event-merging adds a secondary event queue to the simulation, such that some of the newly-arriving events are buffered temporarily in this queue. Whenever a new event arrives at the simulation, the simulation inserts it in one 44 of its event queues, depending on the event-destination-node time bucket. Every newly- arriving event will be queued in the secondary queue unless its destination-node time bucket is zero. In that case, the event is queued in the primary queue directly since this event is not going to be merged. The algorithm carries out the following: 1. Scan the secondary queue events and group these events by their destination nodes. 2. Scan each group of events to determine the sets of events that fall within the par- ticular destination node’s time bucket. 3. Create a meta event that represent each set of events. 4. Queue the meta event in the primary queue, and remove the merged sets of original events from the simulation Figure 3.9 represents the event-merging acceleration algorithm flowchart, and Algo- rithm 4 briefly represents the event-merging algorithm design. The speed up and error analysis for event merging is similar to lookahead event-driven simulation acceleration, and will be performed as part of the future research. Algorithm 4 Event-Driven Simulation with Event Merging 1: Load input events into the queue 2: initialize simulation clock to zero 3: while event queue is not empty do 4: define current event as the event with the earliest delivery time 5: define current neuron as the current event destination 6: carry out a transfer function computation to compute the new state 7: generate new events to notify post-synaptic neurons about the new state 8: add new events to the secondary queue 9: merge events in the secondary queue 10: add merged events to the main queue 11: set simulation clock to the current event delivery time 12: end while 45 The Event-Merging Time Bucket The event-merging time bucket analysis is similar to the lookahead event-driven algo- rithm. However, the event-merging time bucket should be smaller than the lookahead time bucket. The event-merging is carried out before queueing and it is possible that the simulator changes the destination node time bucket to a smaller value later in the simula- tion time before executing the meta event that uses the previous bigger time-bucket value. As such, choosing the current time bucket for event merging may result in bigger accu- racy loss. Thus, choosing a small time bucket by having a small k factor relative to the lookahead algorithm in the time bucket equation is recommended. The event-merging and lookahead simulation algorithms complement each other with the former working before event insertion and the later after event fetching to reduce the queue complexity, and reduce computation rate. 3.4.3 Self-Regulating Event-Driven Simulation Algorithm Overview The premise of the self-regulating event-driven simulation algorithm is to completely eliminate newly arriving low-value events from the simulation to avoid queueing and processing these events. A low-value event is an event that does not change the simulation or it changes it very slightly so that it is not worth consideration. Given a simulation of two nodes A and B, with node A sending events to node B, these are few examples of low-value events: 46 1. The event e i message 3 is exactly equal to the e i1 message. Node B can simply buffer the last value from node A, and assume that this is node A’s current state until it receives an event that invalidates the previous state. 2. The event e i message is almost equal to the e i1 message. Almost equality can mean that messages are within10% of each other, for example. Detecting toler- able percentage depends on the responses of the simulation nodes. The algorithm will trace the receiving nodes’ responses, and measure their sensitivities to the inputs. 3. The event e i message will not change the receiving node’s response. Operation The self-regulating event-driven simulation algorithm operation relies on determining the impact of an event on the simulation. We propose using a sensitivity table for each node to track the response change for every input. Whenever a noden i receives an input, the difference between the new and previous state is computed. The difference value is stored for the recent M events arriving from this particular input to noden i . The value of M will be chosen in a way to determine the impact of an input on the noden i accurately. After recording the influence of each input on noden i , the simulation can make decisions about the influence of events arriving from each input. For example, if events arriving from a certain input changed the node response within1% for the last M events, the simulation should pay less attention to this particular input and attempt to process its events at a lower frequency. Figure 3.10 represents the event-dropping acceleration flowchart. 3 Message in a neural simulator is any change in the neuron’s membrane potential. 47 We will use statistical analysis to determine the speed up and error tradeoffs for vari- ous input profiles 3.4.4 Differences and Similarities Among the Lookahead, Event Merging and Self-Regulating Algorithms These three algorithms attempt to create periods of inactivity in the analog simulation to create an efficient event-driven simulation. An analog system does not have periods of inactivity, as the inputs and outputs changing continuously that could result in massive event-arrival rate that could limit the simulation scalability. However, these three algo- rithms analyze the responses and make decisions to speed up the simulation by merging events, dropping events and executing events in batches. These operations imply creating periods of inactivity, allowing the simulation to move simulation time using dynamic time steps for each single node instead of fixed, small time steps as in time-driven simulators. The lookahead and event-merging are quite similar in the sense that a set of events that falls within a certain period of time is being approximated or merged together. However, there are two major differences between the two algorithms. Firstly, the event-merging time bucket should be smaller than the lookahead event-driven time bucket. Secondly, event-merging is done before queueing events while lookahead is done after queuing events when a node is ready to process its next event. Event-merging helps to reduce the event queue complexity as it operates before queueing. The self-regulating event-driven approach attempts to remove events entirely from the simulation by determining their impact on the recent simulation history. 48 3.5 Established Simulation Acceleration 3.5.1 Event Data Structure Simulating an event’s delay requires storing the event in a data structure and executing it later when it is scheduled to occur. Such data structures are typically implemented as queues. Given the large event arrival rate, we need a very efficient data structure that allows inserting and searching for events as fast as possible, as Rubtina needs to scan this data structure continuously. The event data structure’s main purpose is to store the neural events between the neurons n i and n j forsd ij simulation seconds such thatsd ij is the delay between the two neurons (8 n i , n j 2 ~ N such that ~ N is all neurons in the retinal system). The event data structure should be optimized for event retrieval. Rubtina implements two approaches for the event data structure using a binary tree and a hash table. The trade-off between the two approaches is implementation simplicity and speed such that a binary tree is simple to implement, while the hash table speeds up the search and insertion operations. The binary tree key is the delivery time for each event. At t 0 the root is chosen to be as close as possible to the average delivery time for all events that are found at t 0 to increase the chances of having a balanced tree. As simulation time succeeds, the tree root node itself gets executed and the root node is now shifted to be the node to the right to avoid searching events that happened in the simulation past. This binary tree data structure achieves search and insertion complexity of O(log 2 (n)) such that n is the number of events. The hash table approach is similar to the calendar queue used by the network simulator NS2 [109]. The hash table keys are increments of the current simulation time, and the values are sorted arrays of events that will be executed on or after the timestamp stored on the associated hash table key. That 49 is, each hash key entry is equal to the summation of the previous key entry andt incr . The value on the hash table keyt 0 +t incr is associated with the array of events that will occur during [t 0 +t incr ;t 0 +2t incr ); these events are sorted whenever a new event is added using quick sort. Att v 2[t 0 +t incr ;t 0 +2t incr ), the first event in the array is compared tot v . If the two timestamps are equal, then the event is executed and removed from the array. The hash table lookups and insertions areO(1). Insertions require resorting a small array of events, and that requires negligible computations. Thet incr can be chosen to be close to the mean time-bucket value allowing fetching an entire time bucket of events with very few lookups. 3.5.2 Macro-models For Complex Nodes A macro-model is an implementation of a node’s transfer function that is faster than the original micro-model and captures the node details to a certain extent. Retinal neurons are very complex, and developing macro-models will increase the simulation scalability [56]. Figure 3.11 illustrates the retinal photoreceptor’s response for various light intensi- ties. In the absence of light, the response is stable at40mv. When the photoreceptor is exposed to a light flash, the photoreceptor hyperpolarizes. That is, the membrane potential becomes more negative because of the reduction of positive ions inside the pho- toreceptor. The hyperpolarization leads to lower Ca2 + concentration, which results in the activation of a regulatory mechanism called light adaptation that regulates the membrane potential by bringing it back to the photoreceptor resting potential. Similar photorecep- tor responses have been recorded by many researchers [5, 3, 4]. An example of a cone macro-model will be described in Section 4.6.1. 50 3.5.3 Distributed Event-Driven Simulation Many distributed and parallel simulations strategies have been proposed [60, 29]. Sev- eral design decisions have to be made to distribute an event-driven simulation on several processors include 1. Event Queue: centralized big queue, or multiple smaller queues on every processor, 2. node assignment into various processors, and 3. synchronization among various processors. Distributed architectures must attempt to reduce IO communications to increase the simulation scalability. The node and queue communication is very frequent. Thus, a centralized queue will imply numerous IO operations; as such, having several queues on every machine might be a better approach. Nodes that communicate with each other fre- quently must be placed in the same processor. Initially, node connectivity can be used to make the node assignment. However, as simulation time progresses, the communication history can be used to reassign nodes into more suitable processors. Regardless of how good the node placement technique is, a synchronization process must be present. The synchronization in the context of event-driven simulation algorithm refers to the commu- nication between events and nodes among various processors. 51 3.6 Simulation Testing Testing has been carried out during the entire development life cycle of Rubtina at the neuron level, layer level and at the whole retinal system level for both the baseline simu- lation and the various accelerated versions of the baseline simulation. Neuron and layer- level testing have been carried out by comparing each neuron output or layer output with its biological counterpart output general characteristics. That is, given our study of the biological retina, we expect certain behavior from the retinal neurons we implemented. For example, we expect the horizontal cells to produce a spatially averaged map of poten- tials. We look at the baseline simulation horizontal cells layer output (which is recorded at the end of every time unit of the simulation), and we analyze it. The lack of spa- tial averaging indicates a simulation error, while having uniform spatial averaging across the entire layer indicates the validity of this layer functionality. The same process is carried out to other retinal cells including photoreceptors, center-surround bipolar cells and ganglion cells. We expect the photoreceptors to carry out photosensitive responses, edge detection from the center-surround bipolar cells and high spiking frequency from the ganglion cells receiving stimuli from bipolar cells. The validation can be carried out by comparing the simulation outputs with actual recordings of the retinal cells [89, 31]. This may require highly detailed simulation mod- els, and this is not the central goal of the research. Our main focus is simulation accel- eration, and as long as the simulation is carried out on plausible biological models that produce fairly valid and consistent output, we should be able to carry out the acceleration and determine the performance tradeoffs. 52 Figure 3.6: Lookahead event-driven algorithm flowchart 53 Figure 3.7: Monotonic response versus non-monotonic (within a time bucket) Figure 3.8: The time bucket size influence on simulation accuracy. 54 Figure 3.9: Event-driven algorithm with event-merging flowchart 55 Figure 3.10: Event-driven algorithm with event-dropping flowchart 56 Figure 3.11: Biological retinal photoreceptor response [76] 57 Chapter 4 Baseline Implementation 4.1 Overview The baseline simulation implementation is a classical event-driven simulation algorithm of the biological retina without any acceleration techniques. The simulator has been developed using the Ruby programming language. The simulator will initialize the neural network given the modeler configurations. Then, the modeler can carry out simulation experiments for various inputs. 4.1.1 Inputs and Outputs The simulator input is a sequence of images (GIF image file). The simulator outputs are the ganglion-cell spike trains, and the outputs of the retinal layers showing the layer state. Ganglion-cell spike trains are stored into files, such that a spike from a ganglion cell GC i is represented by the GC i location and spike time (GC ix , GC iy , Time). The outputs of inner layers including the layers of photoreceptor cells, horizontal cells, and bipolar cells are stored into binary images. The inner layers are stored into images whenever the outputs of the majority of a layer of cells change, producing series of images representing the formation of these layers’ potentials. Storing inner layers’ outputs is important for a retinal simulator to allow modelers to use it for various experiments. A modeler interested in studying the visual cortex will need the ganglion cells spike trains for certain stimuli, while retinal prosthetic modelers 58 who are interested in determining the number of cones to allow for facial recognition will be interested in the outer retinal layers’ outputs. 4.1.2 Event Creation From the Input An event represents any message to be communicated among neurons. Whenever a neuron i receives an event, the neuron i state is recomputed. If the neuron i state has changed, a set of events will be created to notify all neuron i ’s post-synaptic neurons about its new state. In Rubtina, the initial events are created from the input. The animated input image is parsed into a set of frames, and a set of events is created for each frame. An event is created for every colored pixel. That is, if the pixel RGB value is not “#FFFFFF”, then an event is created for it since this pixel is part of an object in the visual scene. The events’ delivery times are determined using the time between the input frames. That is, if the input is a two-frame animation with 1:5 seconds between those two frames, then, the first frame’s events will be delivered att = 0, while the second frame’s events will be delivered att = 1:5. These events will be stored into the simulation event queue, and will get delivered to the retinal photoreceptors chronologically. 4.1.3 Cellular-Level Simulation Rubtina is a cellular-level simulation that offers flexible design. Rubtina maintains an object for each neuron in the simulator that receives its inputs and produces its outputs autonomously. Rubtina is more flexible, compared to more “rigid” designs that abstract a large set of neurons as a sub-system, where the modeler cannot change any cell behavior (connections, delay and transfer function). However, the rigid design is usually faster 59 and easier to scale. An example of a rigid design retinal simulation is Virtual Retina. To change the Virtual Retina biological retinal model, the modeler needs to reformulate the mathematical equations, and making significant code changes. 4.1.4 Ruby Programming Language Ruby is a high-level and object-oriented programming language that offers many libraries. Other languages like Java or C++ can be used, but our knowledge of the Ruby language and the passion for the Ruby compact syntax encouraged us to use Ruby. We call our simulator “Rubtina” as a play on the words ruby and retina. We have used many available ruby libraries (“gems”) for image processing, and configuration parsing. 4.2 Configuration One of Rubtina design goals is to allow modelers to create simulations without the hassle of frequent writing code for most of the time or learning a new simulator syntax to make creating simulation experiments easier and faster. Rubtina is capable of creating simula- tions from the modeler-defined variables in the configuration files and allow modelers to input any specification they desire about their retinal models. Rubtina is ready to accept various retinal attributes that could be desired by modelers. We attempted to make the simulator as configurable as possible, but in certain scenar- ios, the modeler may need to add classes to carry out certain advanced experiments. For example, changing the transfer function for a certain retinal cell or adding a new retinal cell will require code development as explained in Section 4.5.2. Simulation platforms vary in the way they allow modelers to create simulations. Some simulators expect the modeler to use known programming languages like C++, 60 or Python to write his simulations [41], while others expect the modeler to use a simu- lator’s specific syntax [100]. Virtual Retina [117] does not offer any biomimetic way to configure its simulation; the modeler is expected to determine various spatial and tempo- ral filters parameters for his retinal model. Rubtina design does not assume any particular behavior of the biological retina that is hardcoded into Rubtina. Rubtina is fully configurable to offer high modeling flexibility and hence we refer to Rubtina as a simulation framework instead of a simulator. We studied the neuroscience literature to determine the attributes that we need to consider to allow creating system-wide simulations on Rubtina. These parameters are discussed in the following sections. 4.2.1 Generic Neural Configurations Most retinal neurons share many attributes. Most of these attributes are configurable without the need to write code, and include the following parameters: 1 - Location (x, y, z): every neuron must have a spatial location in the retinal system. The retinal system is a 3D system since the retinal layers are stacked above each other. Thus, a spatial location using the three-coordinate system is necessary. 2 - Transfer function: whenever a neuron receives an event, a new state computation is carried out using the modeler-specified transfer function for that neuron. Rubtina provides an implementation for the major retinal neurons, and allows a methodical approach to add new neural implementations as explained in Section 4.5.2. 3 - Connectivity: a representation of the synaptic connections of the neuron with other neurons. This entails various configuration parameters including: 61 1 - post-synaptic neurons, 2 - connectivity patterns: Retinal neurons’ synaptic connections vary in terms of direction and length across the retinal system. In terms of direction, the retina has feedback, feed forward or lateral synaptic connections. In terms of length, the retina has short synaptic connections with surrounding neurons (horizontal cells), and long connections across the retinal system (wide-field amacrine cells), 3 - synaptic connection type: inhibitory or excitatory, and 4 - synaptic connection delay. 4.2.2 Retinal Layer Configurations Retinal neurons are stored in a set of layers. Each layer is implemented using a two- dimensional array of neurons. The neuron location inside a particular layer is represented by its two-dimensional location (x, y) and the layer location on the z-axis. Rubtina’s layer concept is a counterpart for the biological retina layer, which contains a set of layers stacked above each other. Layers configurations include: 1 - Region: an area of interest within a layer that can be used to set certain attributes for the region neurons. Retinal neurons vary across the associated layer in different aspects including their density. For example, the retinal cones’ and rods’ densi- ties vary exponentially across the photoreceptor layer as the distance increases or decreases from the center of the fovea. 2 - Neurons types, and locations and their densities in the layer or its regions. 62 4.2.3 Neuron Specific Configurations Rubtina offers several built-in neurons for the major retinal neurons. These neurons include the configuration listed in Section 4.2.1. However, each neuron behavior differs than other neurons demanding additional specific configuration for each neuron, includ- ing the following: 1 - Photoreceptors: a. Cones’ and rods’ percentages of the total photoreceptor count, b. Red, Green and Blue cones’ percentages of the total cone count, c. Density across the photoreceptor layer represented by a set of regions (hollow circles), d. Spectral sensitivity, and e. Width of the photoreceptor outer segment. 2 - Horizontal cells: modeler can configure gap-junction behavior with neighboring cells including: a. Gap junction reach: maximum distance at which two horizontal cells can form a gap junction, and b. Percentage of the cell potential communicated with neighboring horizontal cells. 3 - Bipolar cells: a. Connectivity patterns to the center including location of pre-synaptic neurons, their spatial distribution (distal or close) and their count, and 63 b. Connectivity patterns to the surround. 4 - Ganglion cells: a. Spiking frequency for each input combination on the center and surround, and b. Connectivity patterns to the center and surround. 4.3 Cells to Layer Assignment 4.3.1 Photoreceptor Layer Cell Assignment The retinal cones and rods are distributed non-uniformly across the photoreceptor layer. The fovea center contains the highest cone density, and lowest rod density. The density of cones and rods increase and decrease respectively proportional to the distance from the fovea center. Figures 4.1 and 4.2 illustrate the cones’ and rods’ densities across the photoreceptor layer given the radius from the fovea’s center [55]. The cones’ and rods’ distributions can be viewed as a set of rings with various densi- ties as in Figure 4.3, such that each ring is characterized by its density, and the radiuses of the outer and inner circles of the ring. To create this non-uniform distribution, let us imagine that we have a two-dimensional matrix of empty “spots” as in Figure 4.4. Each spot represents a possible photoreceptor location. The matrix spot density is equal to the maximum density of the cells in the photoreceptor layer. For example, the cone’s matrix will have the fovea’s center cone density as the spot density (125k/mm 2 for the human retina [55]). Each spot may or may not have a photoreceptor cell inside it. The cell assignment into these spots depends on the each cell’s biological density in the ring that contains the spot. For example, the fovea center cone’s density is 125k/mm 2 , 64 Figure 4.1: The cones’ density given the radius from the fovea’s center [55]. which is the same density as the matrix. Then, every spot that falls within the radius of the foveal center contains a cone as illustrated in Figure 4.5. To carry out the cell assignment for the outer photoreceptor density rings, let us pick the ring 0.8mm-1mm. In that ring, the cone density is 10k cells/mm 2 , and the rod density is 18k cells/mm 2 . Using this biological data, we can compute the following: Probability for a spot to have a cone=10/125k=8% Probability for a spot to have a rod=18/125k=14% If R is a random number belongs to [0, 99], then: – 0<R<= 8 the spot contains a cone – 86<R<= 100 the spot contains a rod – 8<R<= 86 the spot is empty 65 Figure 4.2: The rod’s density given the radius from the fovea’s center [55]. Figure 4.3: Cones’ and rods’ distributions [44] Figure 4.6 illustrates the cone density for several rings such that the darker color repre- sents higher cone density. Using a similar approach, we can decide if the cone is Red or Green or Blue. The cones’ percentages in the human’s retina are 75% Red, 20% Green, and 5% Blue [90]. Given that, we can use a random integer variable that ranges between 0 and 99 to deter- mine the cone’s type: 66 Figure 4.4: Empty spot matrix 0<R<75 the cone is Red 75<=R<95 the cone is Green 95<= R<100 the cone is Blue The human’s retina area is approximately around 1100mm 2 . Given our two- dimensional spot matrix with a density of 125k spots/mm 2 , we have a matrix with a 137.5 million spots. The majority of the spots in this matrix are empty. To avoid creating this matrix, we use a sparse matrix. The sparse matrix is implemented using a hash table with the row and column used as a key, and the cell type is the value of that key. If a location represented by a location (x, y) does not have an entry in the hash table, then this location does not have a cell. Figure 4.7 illustrates the output of the photoreceptor layer given the existence of three density rings. 67 Figure 4.5: The fovea’s center cone distribution An alternative approach to map the photoreceptor cells within their layer is to use log polar mapping as illustrated in Figure 4.8 [97]. 4.3.2 Subsequent Layer Cell Assignment Once the photoreceptor-layer cells have been assigned to their locations, the simulator can carry out the cell assignment relative to the assigned photoreceptor cells for the remaining layers using the modeler’s configuration data including connectivity patters and connections statistics. For example, given the modelers’ configurations for the hor- izontal cells including the average number of connections with cones and rods, we can place an appropriate number of horizontal cells in their layer relative to those cones and rods. Once the horizontal cells are assigned, we can carry out the bipolar cells layer assignment relative to both horizontal cells and photoreceptors, and so on. 68 Figure 4.6: Cone density of several density rings Figure 4.7: Photoreceptor layer output given three rings 4.4 Baseline Operation Given a certain retinal configuration, the baseline simulation will go through the follow- ing steps: 1 - Simulation initialization: 69 Figure 4.8: Log polar mapping of the photoreceptor cells 1 - Configuration parsing. 2 - Neuron instantiation: using configuration data, neurons are created and stored in various data structures. 3 - Neuron placement: a location is assigned to each neuron, and its dendrites. 4 - Connectivity: neurons are linked together based on the modeler’s configura- tion. 5 - Event queue initialization. 6 - The input movie is parsed into frames, and a set of events is created to repre- sent each frame. These events will be added to the event queue. 2 - Run time: 1 - The simulator searches the event queue for the event (e i ) with the earliest delivery time. 2 - For each earliest event e i : 1 - The current simulation time is set to be the e i delivery time. 70 2 - The event e i destination neuron n i is determined. 3 - Neuron n i processes the incoming event e i : a new state computation is carried out for neuron n i . 4 - Neuron n i generates new events: if the neuron n i state has changed, then a set of events will be generated to communicate the new state to post- synaptic neurons. The new event’s delivery time is equal summation of current simulation time and n i ’s synaptic delay. 5 - Newly generated events are added the queue. 6 - Event e i is removed from the event queue. 3 - The simulation will end when the event queue is empty. 3 - Output analysis 4.5 Software Architecture 4.5.1 Event-Driven Simulation Strategy In an event-driven simulation strategy, each neuron communicates its state to post- synaptic neurons using events. These events are stored in a temporary data structure called the event queue. The simulator fetches the event with earliest delivery time, and processes it. New events will get produced during simulation, and will be added to the event queue. We discussed the reasons to select an event-driven simulation strategy over a time-driven strategy in Section 3.2. Algorithm 5 illustrates the baseline event driven simulation algorithm. 71 Algorithm 5 The Baseline Event-Driven Simulation Algorithm 1: queue=input 2: sim clock =0 3: while queue is not empty do 4: event=queue:readfread the event with the earliest delivery timeg 5: current neuron = event:destination neuronfcurrent neuron is the current event destination neurong 6: queue:delete(event)fdelete current event from the queueg 7: current neuron:process incoming event(event)fcarry out a new state compu- tation for current neuron given the arriving eventg 8: new events = current neuron:generate output eventsfgenerate new events to notify post-synaptic neurons about the new stateg 9: queue=queue+new eventsfadd the new events to the queueg 10: sim clock =event:delivery timefsimulation clock is set to the delivery time of the current eventg 11: end while 4.5.2 Neuron Base Class and Extensibility All baseline retinal neurons must inherit a parent “Neuron” class that defines the Rubtina neural class structure. For example, a retinal neuron without the methods process incoming event and generate output events will not function at all in the baseline simulation since the entire simulation assumes their existence. Employing object-oriented concepts (polymorphism and inheritance) in the baseline simulation allows the modeler to extend the Rubtina neural library easily since the mod- eler does not need to edit the simulation code, but he needs to inherit the Neuron class and develop the new cells given Rubtina specifications. Algorithm 5 utilizes polymorphism, as an instance of the generic Neuron class is used to call the “child” class function to process events, and generate new events. The polymorphism requires every child class to inherit the parent Neuron class. In Algorithm 5, the simulation determines the current event destination neuron, then calls the following methods from the neuron object: 72 1 - process incoming event: processes the current event and will result usually in an update for the neuron state, 2 - generate output events: generates new events if the current neuron state has changed, and the new events will be stored in the event queue using their deliv- ery time as keys. Since every cell in the baseline simulation inherits the Neuron base class, an object from the type Neuron will allow calling the child cell methods. 4.5.3 Event Queue Both the input events and newly generated events during simulation are stored in the event queue data structure. These events will be fetched chronologically during the sim- ulation to be processed. Different data structures can be used including a linked list, a binary tree and a hash table. The trade-offs of each data structure include implementa- tion complexity, used memory and lookup speed (find the event with the earliest delivery time). The baseline default event queue is a calendar queue built using a hash table, such that the table keys are delivery times and their associated values are arrays of event objects to be processed at that particular delivery time as illustrated in Figure 4.9. In Section 3.5.1, we discussed two implementations and their trade offs. 73 Figure 4.9: Calendar queue 4.6 Built-in Neural Models 4.6.1 Photoreceptor Cells - Macro-model A macro-model implementation is a simple implementation of the system to increase the system scalability without significant accuracy loss. Many macro-models have been pro- posed [69, 93]. The cone macro-model is basically a mathematical equation that approx- imates the response of the photoreceptor given various parameters including input light intensity and photoreceptor outer segment width. Photoreceptors responses can be approximated using a simple formula that uses two linear functions. The first linear function models the photoreceptor transduction, and the second function models the light adaptation. Both functions are function of light intensity and time. We wrote a simple Java program to render these linear functions for various light intensities. We use the grayscale value as an indication of the light intensity. Figure 4.10 illustrates the program output for four intensities (30, 100, 150 and 255). The cone model used can be described using the following mathematical function: f(t;intensity)= 8 > < > : ((0:42(intensity+1))=255)t40 tt adaptation (0:03t+(10=intensity)65 t>t adaptation We modeled the responses in Figure 6.5 as a function that consists of two linear 74 Figure 4.10: A simple cone macro-model functions for the periods before and after the adaptation time 1 . We chose one of the responses to compute the linear function’s parameters, and then modified the linear function to accommodate for intensity such that each function produces reasonable results for various intensities. Notice that for the hyperpolarization linear function, the intensity was multiplied by the time since it changes the slope of the response. For the light adaptation function, the light intensity was shifting the response to a higher or lower level and this is why it was added to the linear function. The neuroscience literature includes more accurate macro-models [69, 93]. The fol- lowing formula represents the cone transduction macro-model: V transduction (t)=V max (1e tK ) Such that: V max : the maximum possible response. t: time K: represents the width of the photoreceptor segment that absorbs photons We carried out Matlab simulations for various K values, and the response looks similar to the biological counterpart as illustrated in Figure 4.11. 1 The adaptation time (t adaptation ) is the time at which the photoreceptor starts depolarizing. t adaptation varies for different intensities and could be described as a function of intensity. 75 Figure 4.11: Cone macro-model adaptation response The following formula is the cone adaptation macro-model [69]: V adaptation (t)=S(t n e t=cone ) n: number of series biological processes that occur during rebuilding the cone photopig- ments. cone : photopigments recovery time constant S: scaling factor The simulation for the adaptation responses is illustrated in Figure 4.12. And, the final cone macro-model equation is: V cone (t)= 8 > < > : V max (1e tK ) tt adaptation S(t n e t=cone ) t>t adaptation 76 Figure 4.12: Cone macro-model transduction response 4.6.2 Horizontal Cells Horizontal cells carry out a continuous spatial averaging process within their cells and through their gap junction connections. Once horizontal cells receive inputs from the photoreceptor, the spatial averaging is carried out, then they inhibit the photoreceptors to contribute to the adaptation process. Essentially, horizontal cells are subtracting the “DC” factor from the input light. Each horizontal cell receives events from photoreceptors and neighboring horizontal cells. Events from horizontal cells will result in carrying out averaging, while photore- ceptor events will change the horizontal cell potential to a value similar or close to their potential. The gap junction communication increases the simulation mean number of events very rapidly due to the continuous state “broadcasting” to neighboring horizontal cells. 77 4.6.3 Center-Surround Bipolar Cells Bipolar cells carry out a straightforward center-surround computation. These cells receive inputs from both the photoreceptors and horizontal cells on their center and sur- round respectively. The center-surround computation will result in contrast enhancement and edge detection. 4.6.4 Spiking Ganglion Cells An ON center-surround ganglion cell will spike very frequently if there is stimulus on the center, and no stimulus on the surround which means that the pixel associated with this particular center has higher light intensity than surrounding pixels [50]. Figure 4.13 illustrates the center-surround ganglion cells behavior and Tables 4.1 and 4.2 illustrate the pixel relative light intensity, given spike frequency for the ON and OFF center-surround ganglion cells respectively. The modeler will determine the spiking frequency of the ganglion cell for both the ON and OFF ganglion cells for the following cases: 1 - Stimulus on both the center and surround 2 - Stimulus on center only 3 - Stimulus on surround only 4 - No stimulus Rubtina will use the spiking frequencies to create the desired spike trains. 78 Figure 4.13: Ganglion cell center-surround behavior. Image is taken from [84] and mod- ified. Spike frequency ON Ganglion Cell None Surround pixels have higher intensity Low Center pixel intensity similar to its surrounding pixels High Center Pixel has higher intensity Table 4.1: ON ganglion cell behavior 4.7 Experiments 4.7.1 The Model We carried out simulations for various inputs for a retinal model illustrated in Figure 4.14. The model includes the major retinal cells. We have carried out simulations in Sections 4.7.2 and 4.7.3 given the following assumptions: Spike frequency OFF Ganglion Cell None Center Pixel has higher intensity Low Center pixel intensity similar to its surrounding pixels High Surround pixels have higher intensity Table 4.2: OFF ganglion cell behavior 79 1 - Input frame size determines photoreceptor count (1 to 1 relationship) 2 - Input frames are loaded into the event queue after creating an event for each colored pixel 3 - The grayscale value of each pixel represents the neuron state 4 - Models used are simple but capture general behavior Figure 4.14: Baseline model (PR: photoreceptor, HC: horizontal cell, BC: bipolar cell, GC: ganglion cell) The delays among various retinal cells have to be chosen carefully to produce valid results. For example, the delay between the cone and bipolar cell must be approximately equal to the delay between the cone and horizontal cell, and the delay between the hor- izontal cell and bipolar cell. Having these values will allow the bipolar cells to carry 80 out the contrast enhancement for the same visual scene “frame” by subtracting its cone outputs from its spatially averaged output (represented by the horizontal cells output). Running the Simulation The simulation is available on GitHub. To run the simulation on any machine, follow these steps: 1 - Install ruby programming language 2 - Install “RMagick” gem: an image processing package [27] 3 - Install Git client: Git is the used source control 4 - Clone Baseline repository using the following command: git clone git@github.com:leadicious/thesis simulation.git 5 - Run the simulation using the following command: ruby the- sis simulation/simulation.rb 6 - Outputs will be stored into the folder thesis simulation/output 4.7.2 Simulation Output for a Square Image Input Figure 4.15 illustrates the outputs of the photoreceptors, horizontal cells and bipolar cells for a static input red square image. No color vision has been implemented, so the pho- toreceptors are showing the same input but in grayscale. Both horizontal cells and bipolar cells carry out their tasks properly as Figure 4.15, such that horizontal cells are perform- ing spatial averaging, and bipolar cells are performing edge detection. The horizontal cells communicate their states to neighboring cells with “weaker” state (lower membrane potential). This results in the expansion of the square area slowly as 81 Figure 4.15: Simulation output for a square potential is being shared from the dark area to the surrounding area in the image as illustrated in Figure 4.16 Figure 4.16: Horizontal cells output expansion Ganglion cells spike trains have been recorded into a file. Their output includes the location of the cell, and the spike time, and this output looks like: [(30;30;2:5);(30;31;2:5);(30;32;2:5);(30;33;2:5);(30;34;2:5);:::] The ganglion cells are center surround structures, and will spike with various frequencies given their center and surround inputs as illustrated in Table 4.2. The output above illustrates that ganglion cells under one of the square input image edges is firing indicating the existence of an “interesting” piece of information. 4.7.3 Simulation Output For a Letter “D” Image Input Figure 4.17 illustrates the simulation output for the capital letter D input for photorecep- tors, horizontal cells and bipolar cells. 82 Figure 4.17: Simulation output for the capital letter D 83 Chapter 5 Experiments 5.1 Overview 5.1.1 Retinal Model The retinal model used is shown in Figure 5.1. The model includes photoreceptors, hori- zontal cells, center-surround bipolar cells and center-surround ganglion cells. The model does not include the amacrine cells. The delay between a photoreceptor and a center- surround horizontal cell is two time units, and the delay between any other two neurons is one time unit. Neural models used are simple but capture the general behavior of each cell. Rubtina is designed and implemented in a way to allow modelers to change the neural model behavior in fairly a simple way as explained in Section 4.5.2. 5.1.2 Horizontal Cells Layer Propagation Horizontal cells continuously communicate their potentials to all neighboring horizontal cells using gap junctions, creating a highly active communication system. To limit the activity in this layer, the simulation offers a way to reduce this activity using the variable “LHCP” (Limit Horizontal Cell Layer Propagation). When LHCP=3, an event produced from a horizontal cell will propagate to neighboring cells for three cells only. The LHCP value for that event will be decremented every time it reaches a destination, and once 84 Figure 5.1: Retinal model LHCP=0, the event is ignored and removed from the event queue. By default, LHCP is disabled (LHCP =1). Enabling LHCP allows modelers to carry out more advanced experiments faster with- out incurring major accuracy loss especially in major output characteristics like edges. Horizontal cell events propagation directly influences the horizontal cells’ spatial averag- ing functionality, which contributes to the contrast enhancement and edge detection func- tionalities by the bipolar cells. Using our simulation experiments, we found that allowing the horizontal cells to propagate three to five times will produce visually acceptably out- put with recognizable edges, and recognizable spatial averaging. 85 5.1.3 Simulation Error Computation Simulation error is computed by comparing the acceleration algorithm results for all lay- ers with the baseline simulation results. That is, the output for each neuron in each layer is compared with the same baseline output. Whenever a difference is found between the accelerated versions and baseline, the pixel is considered “erroneous”. If the difference value is “high” (e.g. larger than 0:9Vmax), the error is considered “highly erroneous”. We compute the percentage of errors relative to the number of outputs, and the number of highly erroneous outputs relative to the errors count. We validate the baseline output by visually inspecting the various retinal layers’ outputs, and checking that these layers are performing their functionalities. For example, the horizontal-cell layer must be perform- ing spatial averaging to be considered valid and the center-surround bipolar-cell layer must be performing edge detection to be considered valid. Validating the model can be carried out by comparing the outputs with the biological counterpart, but our goal is to create a fairly functional retinal model that allows us to test our simulation acceleration algorithms, and the biological accuracy of our approach should suffice for our purposes. 5.1.4 Acceleration Algorithm Computation Cost The self-regulating acceleration algorithm computation cost is negligible since the algo- rithm is simply ignoring creating events when the state-change percentage is smaller than a certain threshold determined by the modeler. The lookahead acceleration algorithm computation cost includes the following: 1 - Fetching events with delivery times in the period [t current , t current +TB] for the cur- rent event assuming the current event’s time bucket is greater than zero. Rubtina’s 86 event queue is a calendar queue. This queue keys are delivery times of the simula- tion events, and the values are “sub-arrays” of events that will be delivered at the delivery time in the associated key. Given this data structure, the lookahead accel- eration algorithm will search several sub-arrays to fetch the time bucket events. Thus, the search complexity is m:O(s), such that m is the number of searched sub-arrays, ands is the length of the sub-array. 2 - Piecewise linear approximation computation: requires two transfer function com- putations of the current event destination neuron. The piecewise linear approxima- tion involves linear computations using the last two computed neural states, and linear computational complexity is fairly negligible. 3 - Storing and searching for approximated responses for the current event destination neuron. The responses are stored in a hash table that offers constant search and store operations. The event-merging acceleration algorithm complexity includes 1 - Fetching time-bucket events: the computational complexity is the same as fetching the time-bucket events for lookahead. 2 - Deleting time-bucket events from the queue and inserting the new meta event. Deleting an event inside a Rubtina calendar queue requires searching for its key (constant complexity) and searching the sub-array for the event (O(s)). Thus, the complexity of deleting one event is O(C)+O(s). Inserting an event requires search- ing for its sub-array (O(C)). 87 5.2 Letter “L” Experiments We carried out several simulation experiments using inputs that contain an image of the letter “L”. The inputs have been selected carefully to exploit the strengths and weaknesses of the proposed acceleration algorithm. 5.2.1 One Frame With Letter “L” Baseline We carried out the baseline simulation for the letter “L” image illustrated in Figure 5.2 disabling and then enabling LHCP. Table 5.1 illustrates the queue lengths at various simu- lation times. When LHCP is disabled, the event queue continues to have events due to the continuous communication at the horizontal cell layer. However, when LHCP is enabled, the horizontal-cell layer stops producing events and the simulation ends at t=6. Figure 5.3 illustrates the baseline event queue growth. Note that the event growth dropped at t = 3 since the center-surround bipolar cells carried out their computation at that time resulting in removing the events received from both photoreceptors and horizontal cells during the previous two time units. Aftert = 3, the events are mainly created due to the horizontal cells continuous communication. Lookahead We carried out the lookahead acceleration algorithm at the horizontal cell layer. Table 5.2 and Figure 5.4 illustrate the event queue growth for baseline and lookahead withTB HC = 2 and TB HC = 4, QL is abbreviation for queue length. The lookahead acceleration algorithm performed exactly like baseline since it needs three or more points (each point 88 Figure 5.2: Image of the capital letter “L” Time Queue length (LHCP = false) Queue Length (LHCP = true) 1 1515 1515 2 2239 2239 3 1492 1492 4 2084 2084 5 2034 2034 6 2125 0 7 2221 0 8 2317 0 9 2413 0 10 2509 0 11 2605 0 Table 5.1: Baseline queue lengths with horizontal-cell propagation limitation turned off and on for the one frame with letter “L” input constitute of a distinct time, and a response value) to carry out a lookahead computation. However, given the uniform one-time-unit delay of the horizontal cell layer, the necessary number of points has never been available and no acceleration has been carried out. Event-Merging We carried out the event-merging acceleration algorithm at the horizontal cell layer. Table 5.3 and Figure 5.5 illustrate the queue lengths of the baseline and the event-merging 89 Figure 5.3: Baseline event growth for the one frame with letter “L” input simulation algorithm. The event queue growth was slower than baseline when the event- merging acceleration algorithm was enabled since it scans the event-queue and merges events with the same destination neuron that falls within a certain time bucket. The merging process reduces the event growth as illustrated in Table 5.3, but changing the exact delivery time of events will produce an error. The achieved speed up is 13% and error is14%. Self-Regulating We carried out the self-regulation algorithm at the horizontal cell layer. In the self- regulating algorithm, incoming events are ignored if their potential is below a certain threshold value determined by the modeler. The term “HC STOP” refers to the small- est event potential that can be accepted by a horizontal cell. When HC STOP is larger, the event queue growth is limited since a relatively high HC STOP value will result in 90 Time Baseline QL Lookahead QL (TB HC =2) Lookahead QL (TB HC =4) 1 1515 1515 1515 2 2239 2239 2239 3 1492 1492 1492 4 2084 2084 2084 5 2034 2034 2034 6 2125 2125 2125 7 2221 2221 2221 8 2317 2317 2317 9 2413 2413 2413 10 2509 2509 2509 11 2605 2605 2605 Table 5.2: Baseline and lookahead acceleration algorithm queue lengths for the one frame with letter “L” input. QL is the queue length Figure 5.4: Baseline and the lookahead acceleration algorithm queue lengths for the one frame with letter “L” input 91 Time Baseline QL Event-Merging QL (TB HC ) = 2 1 1515 1515 2 2239 2239 3 1492 1215 4 2084 1766 5 2034 1617 6 2125 1919 7 2221 1790 8 2317 2083 9 2413 1966 10 2509 2247 11 2605 2142 Table 5.3: Baseline and the event-merging acceleration algorithm queue lengths for the one frame with letter “L” input Time Baseline QL QL (HC STOP = 0.90Vmax) QL (HC STOP = 0.88Vmax) 1 1515 1515 1515 2 2239 2239 2239 3 1492 920 1255 4 2084 1151 1715 5 2034 1293 1794 6 0 0 0 Table 5.4: Baseline and self-regulating algorithm queue lengths for several HC STOP values for the one frame with letter “L” input dropping more events. Event queue lengths are illustrated in Table 5.4 and Figure 5.6. Error and speed up computations are illustrated in Tables 5.5 and 5.6 and Figure 5.7. The algorithm performed well for this input by increasing the simulation speed faster than the error increase. LHCP has been enabled for these experiments. When HC STOP value is larger, more events are ignored and the queue length at various times is smaller as illustrated in Table 5.4, however, ignoring events with larger HC STOP value produced higher error as illustrated in Table 5.5. 92 Figure 5.5: Baseline and the event-merging acceleration algorithm queue lengths for the one frame with letter “L” input HC STOP Outputs Erroneous Outputs Highly Erroneous % Mean Error Error % 0.90Vmax 85406 8689 98% 0.31Vmax 10% 0.88Vmax 85406 2520 95% 0.36Vmax 2% Table 5.5: Self-regulating algorithm error for the one frame with letter “L” input HC STOP Total Event Count Speed Up Error 0 (Baseline) 9364 0% 0% 0.90Vmax 7118 23.9% 10% 0.88Vmax 8518 9% 2% Table 5.6: Self-regulating algorithm performance for the one frame with letter “L” input 93 Figure 5.6: Baseline and self-regulating algorithm queue lengths for several HC STOP values for the one frame with letter “L” input 5.2.2 Ten Frames With Letter “L” With Linear Intensity Increase Given a movie input as illustrated in Figure 5.8 that consists of ten frames, the intensity of each pixel is represented by the function: Intensity(t)=18t+20. The time between any two frames is one time unit. Lookahead The photoreceptor model response is directly proportional to the input intensity. Given that the input is linearly changing and the lookahead acceleration approximates the node’s future responses linearly, we carried out the lookahead acceleration algorithm at the pho- toreceptor layer. Event queue growth is illustrated in Table 5.7 and Figure 5.9. Unlike the 94 Figure 5.7: Self-regulating algorithm performance for the one frame with letter “L” input Figure 5.8: Sequence of letter L images with linear intensity increase (the actual input contains one input image for each frame) event-merging and the self-regulating algorithms, the event queue length is the same for the lookahead experiments, and this is ideally what is expected from the lookahead. The lookahead acceleration algorithm approximates event responses for a neuron using its last two responses. The approximated responses are stored in a hash table keyed by delivery times. Once an event arrives to this neuron with an approximated response stored in the hash table, the algorithm uses the stored approximated response value from the hash 95 Time Baseline QL Lookahead QL (TB PR =10) Lookahead QL (TB PR =5) 1 6060 6060 6060 2 3249 3249 3249 3 3778 3778 3778 4 3706 3706 3706 5 4440 4440 4440 6 1248 1248 1248 7 1010 1010 1010 8 1010 1010 1010 9 1010 1010 1010 10 1010 1010 1010 11 1010 1010 1010 12 0 0 0 Table 5.7: Lookahead queue lengths for the letter “L” movie with linear intensity increase input table and removes the event from the event queue without carrying out a transfer func- tion computation. This process should not involve any queue length modification unless the approximations are not accurate. Table 5.8 illustrates the error and speed up. Table 5.9 and Figure 5.10 illustrates the lookahead acceleration algorithm performance for this input in terms of the number of transfer function computations and approximations. The lookahead acceleration algorithm performed very well for this particular example since the input is linearly increasing and we are performing the lookahead algorithm on the photoreceptor layer. The photoreceptor linear response is also the main reason for having the same speed up and error numbers for both time-bucket values since the approximation of a linear response for the neuron next 10 seconds will produce the same results whether the approximation is carried out once (TB PR =10) or twice (TB PR =5). LHCP has been enabled for these experiments. The lookahead algorithm speed up is the ratio of number of approximations to the baseline number of computations. 96 Figure 5.9: Lookahead acceleration algorithm for the letter “L” movie with linear inten- sity increase input TB PR Outputs Erroneous Outputs Highly Erroneous % Mean Error Error % 5 573004 0 0% 0.00Vmax 0% 10 573004 0 0% 0.00Vmax 0% Table 5.8: Lookahead algorithm error for the letter “L” movie with linear intensity increase input Event-Merging We carried out the event-merging acceleration algorithm on the horizontal cell layer. The delays of the horizontal cells are not uniform in this experiment. Each horizontal cell has connections with delays of 1, 2, and 3 time units. The variable delay will illustrate the event-merging algorithm’s time-bucket value influence on speed up and error trade- offs. Table 5.10 and Figure 5.11 illustrate the event queue growth, and Table 5.11 and 97 Time Bucket Computations Approximations Speed Up Error 0 (Baseline) 26521 0 0% 0% TB PR =5 21471 5050 19.04% 0% TB PR =10 21471 5050 19.04% 0% Table 5.9: The lookahead acceleration algorithm performance for the letter “L” movie with linear intensity increase input Figure 5.10: The lookahead acceleration algorithm performance for the letter “L” movie with linear intensity increase input Figure 5.12 illustrate the performance of the algorithm. Notice that whenTB HC =4 the algorithm merged 5% more events than when TB HC =2 since with larger time bucket, the algorithm fetched more events and brought it earlier. The acceleration algorithm increased the performance twice as fast as the error increase. 98 Time Baseline QL Event-Merging QL (TB HC = 2) Event-Merging QL (TB HC = 4) 1 6060 6060 6060 2 3249 3249 3249 3 1657 1909 2101 4 2236 2542 2823 5 2391 1870 1970 6 3198 2896 3444 7 3356 2160 2401 8 4112 2792 3573 9 4111 2630 2848 10 5741 3033 3488 11 4492 2363 2527 12 4159 2429 2079 13 564 1181 988 14 1493 1236 1055 15 507 702 543 16 1487 1190 1114 Table 5.10: Event-merging queue lengths for the letter “L” movie with linear intensity increase input 99 Figure 5.11: Event-merging queue lengths for the letter “L” movie with linear intensity increase input 100 TB HC Outputs Erroneous Outputs Highly Erroneous % Mean Error Error Speed Up % 2 748408 187668 86% 0.39Vmax 25% 49.7% 4 748408 204719 88% 0.39Vmax 27% 50.2% Table 5.11: Event-merging performance for the letter “L” movie with linear intensity increase input 101 Figure 5.12: Event-merging performance for the letter “L” movie with linear intensity increase input Self-Regulating Event queue growth is illustrated in Table 5.12 and Figure 5.13. The algorithm perfor- mance is illustrated in Table 5.13 and Figure 5.14. LHCP has been enabled for these experiments. When HC STOP=0.88Vmax, the simulation speed and error increased at the same rate. However, when HC STOP¿0.88Vmax, the error increased at faster rate than the speed, thus, we should avoid running the simulation with those values or choose another acceleration algorithm for this particular layer. 102 Time Baseline QL SR QL (0.90Vmax) SR QL (0.894Vmax) SR QL (0.88Vmax) 1 6060 6060 6060 6060 2 3249 3249 3249 3249 3 3778 3420 3420 3778 4 3706 3243 3255 3704 5 4440 3643 3801 4316 6 1248 1125 1229 1205 7 1010 1010 1010 1010 8 1010 1010 1010 1010 9 1010 1010 1010 1010 10 1010 1010 1010 1010 11 0 0 0 0 Table 5.12: Self-regulating queue lengths at various HC STOP values (SR: self regulating algorithm) for the letter “L” movie with linear intensity increase input Figure 5.13: Self-regulating queue lengths at various HC STOP values for the letter “L” movie with linear intensity increase input 103 HC STOP Outputs Erroneous Outputs Highly Erroneous % Mean Error Error Speed Up % 0.90Vmax 573004 82551 94% 0.25Vmax 14.4% 6.56% 0.894Vmax 573004 33725 88% 0.42Vmax 5.8% 5.5% 0.88Vmax 573004 3601 66% 0.34Vmax 0.6% 0.6% Table 5.13: Self-regulating algorithm error and speed up for the letter “L” movie with linear intensity increase input 104 Figure 5.14: Self-regulating algorithm error and speed up for the letter “L” movie with linear intensity increase input 5.2.3 Ten Frames With Letter “L” With Non-Linear Intensity Increase Given a movie input that consists of ten frames, and the intensity of each input is non- linearly changing, represented by the following function: Intensity(t)=2t 2 +2. The time between any two frames is one time unit. 105 Lookahead The event queue growth is illustrated in Tables 5.14 and Figure 5.15. The algorithm performance is illustrated in Tables 5.15 and 5.16 and Figures 5.16 and 5.17. LHCP has been enabled for these experiments. When TB PR =2, the lookahead performed just like baseline. With frames coming at a rate of 1 frame/second and TB PR =2, the algorithm will never find a set of events that consists of three or more events to carry out a lookahead. The lookahead searches for events that fall in the period [t first event delivery time ,t first event delivery time +TB PR ]. When TB PR =2, finding the necessary set of events will never happen. Table 5.14 illustrates queue lengths for the baseline simulation and the lookahead sim- ulation at various time-bucket values. The queue length values differed at various times (e.g t=5) indicating possible approximation inaccuracy. When the lookahead approxima- tion produces stronger or weaker response for a certain neuron, then the post-synaptic neurons responses will change, leading to the queue length discrepancies. The lookahead acceleration algorithm performed poorly for the non-linearly changing input as illustrated in Table 5.16. The lookahead performed linear approximation for a set of neurons responding non-linearly, producing larger error than speed up. Event-Merging Event queue growth is illustrated in Table 5.17 and Figure 5.19. Algorithm performance is illustrated in Table 5.18 and Figure 5.18. Variable horizontal cell delay is used in this experiment. Notice that the event queue at certain times was larger when the event- merging algorithm was enabled (e.g. at t=4). It appears that the algorithm produced large 106 Time Baseline QL QL (TB PR =10) QL (TB PR =5) QL (TB PR = 3) QL (TB PR =2) 1 6060 6060 6060 6060 6060 2 3249 3249 3249 3249 3249 3 2940 2940 2940 2940 2940 4 3453 3453 3453 3453 3453 5 3044 3044 3182 3182 3182 6 1035 1035 1173 1173 1035 7 1010 1010 1148 1148 1010 8 1010 1010 1148 1010 1010 9 1010 1010 1148 1010 1010 10 1010 1010 1148 1010 1010 11 0 0 138 0 0 Table 5.14: Lookahead queue lengths for various photoreceptor time-bucket values for the letter “L” movie with non-linear intensity increase input TB PR Outputs Erroneous Outputs Highly Erroneous % Mean Error Error 2 573061 0 0% 0.0 0% 3 573061 62127 16% 0.08Vmax 10% 5 573061 221213 32% 0.10Vmax 38% 10 573061 317326 91% 0.45Vmax 55% Table 5.15: The lookahead acceleration algorithm error for the letter “L” movie with non-linear intensity increase input number of meta events that entailed a lot of events that would have arrived in the near simulation future. 107 Figure 5.15: Lookahead queue lengths for various photoreceptor time-bucket values for the letter “L” movie with non-linear intensity increase input TB PR Computations Approximations Speed Up Error 0 (Baseline) 23821 0 0% 0% TB PR =10 19737 5050 20.4% 55% TB PR =5 19185 5050 20.8% 38% TB PR =3 19919 4040 16.8% 10% TB PR =2 23821 0 - - Table 5.16: The lookahead acceleration algorithm performance for the letter “L” movie with non-linear intensity increase input 108 Figure 5.16: The lookahead acceleration algorithm error for the letter “L” movie with non-linear intensity increase input 109 Figure 5.17: The lookahead acceleration algorithm performance for the letter “L” movie with non-linear intensity increase input 110 Time Baseline QL Event-Merging QL (TB HC =2) Event-Merging QL (TB HC =5) 1 6060 6060 6060 2 3249 3249 3249 3 1448 1651 1657 4 2315 3685 4208 5 1440 1446 1442 6 2475 2056 2251 7 1595 1652 1466 8 2456 2112 1950 9 1535 1602 1500 10 2695 2331 2068 11 1012 769 526 12 1974 1751 995 13 545 649 576 14 1571 1346 1129 15 633 738 602 16 1679 1262 1144 Table 5.17: The event-merging acceleration algorithm queue lengths for the letter “L” movie with non-linear intensity increase input 111 Figure 5.18: The event-merging acceleration algorithm queue lengths for the letter “L” movie with non-linear intensity increase input 112 TB HC Outputs Erroneous Outputs Highly Erroneous % Mean Error Error Speed Up 2 749249 165182 81% 0.33Vmax 22% 29% 5 749249 173788 77% 0.33Vmax 23% 32.2% Table 5.18: The event-merging acceleration algorithm performance for the letter “L” movie with non-linear intensity increase input 113 Figure 5.19: The event-merging acceleration algorithm performance for the letter “L” movie with non-linear intensity increase input Self-Regulating Event queue growth is illustrated in Table 5.19 and Figure 5.20. Algorithm performance is illustrated in Table 5.20 and Figure 5.21. 114 Time Baseline QL SR QL (HC STOP=0.92Vmax) SR QL (HC STOP=0.90Vmax) 1 6060 6060 6060 2 3249 3249 3249 3 2940 1844 2368 4 3453 1919 2497 5 3044 1999 2468 6 3140 1811 2466 7 3236 1903 2649 8 3332 2015 2885 9 3428 2123 3072 10 3524 2230 3261 11 2610 1328 2437 12 2706 1436 2622 13 2802 1544 2814 14 2898 1652 3006 15 2994 1760 3198 Table 5.19: The self-regulating algorithm performance (SR: self-regulating algorithm) for the letter “L” movie with non-linear intensity increase input 115 Figure 5.20: The self-regulating algorithm performance for the letter “L” movie with non-linear intensity increase input 116 HC STOP Outputs Erroneous Outputs Highly Erroneous % Mean Error Error Speed Up % 0.92Vmax 1040089 441279 96% 0.32Vmax 42% 33.5% 0.90Vmax 1040089 401151 95% 0.22Vmax 38% 8.83% Table 5.20: The self-regulating algorithm performance for the letter “L” movie with non-linear intensity increase input 117 Figure 5.21: The self-regulating algorithm performance for the letter “L” movie with non-linear intensity increase input 118 5.2.4 Ten Frames Moving Letter “L” With Linear Intensity Increase Given a 10-frame input, each frame contains the image of the capital letter L. The center of the letter L is moving horizontally to the right every time unit. We carried out our experiment using two variations of this input given the speed of the movement: 1 - c x (t) = c x (0) + t: the center of the letter L moves horizontally to the right 1 pixel/second. Given that the frame arrival rate is 1 frame/sec, the entire input image moves one pixel to the right every second. 2 - c x (t) = c x (0) + 3t: the center of the letter L moves horizontally to the right 3 pixels/second. Lookahead For the first input (c x (t) = c x (0) +t), event queue growth is illustrated in Table 5.21 and Figure 5.22. Algorithm performance is illustrated in Table 5.22 and Figure 5.24. For the second input (c x (t) = c x (0) + 3t), event queue growth is illustrated in Table 5.23 and Figure 5.25. Algorithm performance is illustrated in Table 5.24 and Figures 5.26 and 5.27. The lookahead acceleration algorithm did not produce any error since the intensity of the input is linearly changing. However, it did not perform as well as the experiment in Section 5.2.2. The reason for the lower performance in this experiment is the input movement, as the lookahead acceleration algorithm is carrying out approximations for certain photoreceptors that will not be stimulated in the near simulation future since the input is moving away from these particular photoreceptors. 119 Time Baseline QL Lookahead QL (TB PR =5) Lookahead QL (TB PR =10) 1 6060 6060 6060 2 3249 3249 3249 3 3107 3107 3107 4 3525 3525 3525 5 4727 4727 4727 6 1619 1619 1619 7 1010 1010 1010 8 1127 1127 1127 9 1798 1798 1798 10 2906 2906 2906 11 2785 2785 2785 12 733 733 733 13 570 570 570 14 560 560 560 Table 5.21: The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) Figure 5.22: The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) 120 TB PR Computations Approximations Error Speed Up 0 (Baseline) 33776 0 0% 0% 5 29274 4502 0% 13.3% 10 28898 4878 0% 14.4% Table 5.22: The lookahead acceleration algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) Figure 5.23: The lookahead acceleration algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) 121 Figure 5.24: The lookahead acceleration algorithm error for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) 122 Time Baseline Lookahead QL (TB PR =5) Lookahead QL (TB PR =10) 1 6060 6060 6060 2 3249 3249 3249 3 3478 3478 3478 4 4021 4021 4021 5 6075 6075 6075 6 4197 4197 4197 7 4203 4203 4203 8 4291 4291 4291 9 4909 4909 4909 10 5100 5100 5100 11 4166 4166 4166 12 1831 1831 1831 13 822 822 822 14 666 666 666 Table 5.23: The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) TB PR Computations Approximations Error Speed Up 5 49306 3762 0% 7.1% 10 49138 3930 0% 7.4% Table 5.24: The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) 123 Figure 5.25: The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) 124 Figure 5.26: The lookahead acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) 125 Figure 5.27: The lookahead acceleration algorithm error for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) 126 Time Baseline QL QL (TB HC =1) QL (TB HC =2) QL (TB HC =4) 1 6060 6060 6060 6060 2 3249 3249 3249 3249 3 1489 1489 1588 1766 4 1944 2026 2363 2546 5 2389 2127 1735 2073 6 3232 2622 2448 2565 7 3451 2422 2328 2147 8 4787 3403 3476 3733 9 5045 3763 2945 3138 10 6822 3395 3548 3491 11 5385 3488 3017 2804 12 3823 2525 2455 2413 13 1127 1847 1667 1742 14 950 832 917 876 15 1077 1085 954 954 16 993 779 749 708 Table 5.25: The event-merging acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) Event-Merging For the first input, event queue growth is illustrated in Table 5.25 and Figure 5.28. Algo- rithm performance is illustrated in Table 5.26 and Figure 5.29. For the second input, event queue growth is illustrated in Table 5.27 and Figure 5.31. Algorithm performance is illustrated in Table 5.28 and Figure 5.31. 127 Figure 5.28: The event-merging acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) 128 TB HC Outputs Erroneous Outputs Highly Erroneous % Mean Error Error Speed Up 1 658701 207613 70% 0.29Vmax 31% 48.4% 2 658701 235466 85% 0.39Vmax 35% 51.7% 4 658701 237109 88% 0.42Vmax 36% 53.23% Table 5.26: The event-merging acceleration algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) 129 Figure 5.29: The event-merging acceleration algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) 130 Time Baseline QL QL (TB HC =1) QL (TB HC =2) QL (TB HC =4) 1 6060 6060 6060 6060 2 3249 3249 3249 3249 3 1533 1533 1901 2099 4 2379 2754 2993 3004 5 3642 3394 3261 3377 6 4323 3582 3632 3653 7 4835 3849 3789 3793 8 5039 4023 3725 3887 9 5162 3934 4012 3849 10 5836 4118 3960 3957 11 4640 3161 3071 2948 12 4092 2333 2101 2117 13 1459 1588 1456 1427 14 1426 1146 1029 1001 15 1194 1112 1166 1156 16 1292 1063 931 887 Table 5.27: The event-merging acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) 131 Figure 5.30: The event-merging acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) 132 TB HC Outputs Erroneous Outputs Highly Erroneous % Mean Error Error Speed Up 1 864835 282493 70% 0.31Vmax 32% 47.8% 2 864835 343797 80% 0.36Vmax 39% 51.1% 4 864835 347866 82% 0.37Vmax 40% 52.3% Table 5.28: The event-merging acceleration algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) 133 Figure 5.31: The event-merging acceleration algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+3t) Self-Regulating Event queue growth is illustrated in Table 5.29 and Figure 5.32. Algorithm performance is illustrated in Table 5.30 and Figure 5.33. 134 Time Baseline QL SR QL (0.92Vmax) SR QL (0.90Vmax) SR QL (0.88Vmax) 1 6060 6060 6060 6060 2 3249 3249 3249 3249 3 3107 2365 2788 3107 4 3525 1720 2709 3523 5 4727 2037 3344 4584 6 5563 2315 4143 5412 7 6308 2656 5230 6197 8 6974 3379 6400 6913 9 7254 4247 6758 7224 10 7671 5259 6899 7540 11 6858 5043 6257 6716 Table 5.29: The self-regulating acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) Figure 5.32: The self-regulating acceleration algorithm queue lengths for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) 135 HC STOP Outputs Erroneous Outputs Highly Erroneous % Mean Error Error Speed Up % 0.92Vmax 360655 132375 93% 0.35Vmax 36% 40% 0.90Vmax 360655 105311 90% 0.33Vmax 29% 12.9% 0.88Vmax 360655 14013 81% 0.47Vmax 3% 1.2% Table 5.30: The self-regulating algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) 136 Figure 5.33: The self-regulating algorithm performance for the moving letter “L” movie with linear intensity increase input (c x (t)=c x (0)+t) 137 Algorithm Speed-Up Average Error Average Lookahead 15.3% 11.4% Event-Merging 42.5% 29% Self-Regulating 14.9% 17.2% Table 5.31: Letter L simulation experiments summary 5.2.5 Letter L Performance Trade-off Summary Table 5.31 illustrates the speed up and error averages for various letter L simulation experiments. Tables 5.32, 5.33 and 5.34 illustrate the speed up and error tradeoffs for the lookahead, event-merging and self-regulating acceleration algorithms respectively. The lookahead algorithm performed badly when the neural responses non-linearly changed, and performed very well for linearly changing responses. The event-merging performed well at the horizontal cell layer, reducing the layer activity. The self-regulating algorithm speed up average is a little less than the produced average error, however, the error pro- duced by the self-regulating algorithm did not change the output response features like edges. The error is computed by comparing the accelerated algorithm output with the baseline at every time unit for every neuron. For self-regulating, when an event reaches a certain threshold value, it will be ignored. However, the error computation will still consider the lack of event-propagation as an error, resulting with an error average higher than the speed up. 138 Lookahead Acceleration Algorithm Performance Input Time Bucket Speed Up Error Static Image TB HC =2 0% 0% TB HC =4 0% 0% Linear Intensity TB PR =5 19.04% 0% TB PR =10 19.04% 0% Non-Linear TB PR =3 16.8% 10% TB PR =5 20.4% 38% TB PR =10 20.4% 55% Moving (c x (t)=c x (0)+t) TB PR =5 13.3% 0% TB PR =10 14.4% 0% Moving (c x (t)=c x (0)+3t) TB PR =5 7.1% 0% TB PR =10 7.4% 0% Table 5.32: The lookahead acceleration algorithm performance for a set of experiments The Event-Merging Acceleration Algorithm Performance Input Time Bucket Speed Up Error Static Image TB HC = 2 13% 14% Linear Intensity TB HC = 2 49.7% 25% TB HC = 4 50.2% 27% Non-Linear TB HC = 2 29% 22% TB HC = 5 32.2% 23% Moving (c x (t)=c x (0)+t) TB HC = 2 51.7% 35% TB HC = 4 53.23% 36% Moving (c x (t)=c x (0)+3t) TB HC =2 51.1% 39% TB HC =4 52.3% 40% Table 5.33: The event-merging acceleration algorithm performance for a set of experi- ments 139 Self-regulating Acceleration Algorithm Performance Input HC STOP Speed Up Error Static Image 0.90Vmax 23.9% 10% 0.88Vmax 9% 2% Linear Intensity 0.90Vmax 14.4% 6.56% 0.89Vmax 5.8% 5.5% 0.88Vmax 0.6% 0.6% Non-Linear 0.92Vmax 33.5% 42% 0.90Vmax 8.83% 38% Moving (c x (t)=c x (0)+t) 0.92Vmax 40% 36% 0.90Vmax 12.9% 29% 0.88Vmax 1.2% 3% Table 5.34: The self-regulating acceleration algorithm performance for a set of experi- ments 140 5.3 Letter “G” Experiments The letter “G” image is a little more advanced than the letter “L” image since it has rounded areas with pixels that have fairly less dark color as explained in Figure 5.34. In the next few sections, we carried out simulations on various inputs that include the letter G with intensities changing linearly or non-linearly. However, this intensity change cannot be applied to the rounded areas. The lookahead acceleration algorithm will carry out its approximation for every input assuming linear intensity change resulting in higher errors due to the input rounded areas. The following examples illustrate the importance of enabling and disabling the simulation algorithm dynamically for each output for certain period of time as explained in Section 6.2.1. Figure 5.34: Letter G rounded area 5.3.1 Lookahead For Ten Frames With Letter “G” With Linear Intensity Increase Event queue growth is illustrated in Table 5.35 and Figure 5.35. Algorithm performance is illustrated in Table 5.36 and Figure 5.36. 141 Time Baseline QL Lookahead QL (TB PR =5) Lookahead QL (TB PR =10) 1 11731 11731 11731 2 6655 6655 6655 3 9683 9683 9683 4 9422 9406 9406 5 10582 10507 10507 6 11398 11337 11337 7 11351 11329 11329 8 11434 11355 11365 9 11891 11982 11884 10 12625 12678 12617 11 11199 11252 11225 12 2990 2991 2990 13 3023 3023 3023 14 3119 3119 3119 15 3215 3215 3215 Table 5.35: The lookahead acceleration algorithm queue lengths for the letter “G” movie with linear intensity increase input Time Bucket Speed Up Error TB PR =5 7.0% 5% TB PR =10 7.7% 10% Table 5.36: The lookahead acceleration algorithm performance for the letter “G” movie with linear intensity increase input 5.3.2 Lookahead For Ten Frames With Letter “G” With Non-Linear Intensity Increase Event queue growth is illustrated in Table 5.37 and Figure 5.37. Algorithm performance is illustrated in Table 5.38 and Figure 5.38. 142 Figure 5.35: The lookahead acceleration algorithm queue lengths for the letter “G” movie with linear intensity increase input 5.3.3 Lookahead For Ten Frames Moving Letter “G” With Linear Intensity Increase Event queue growth is illustrated in Table 5.39 and Figure 5.39. Algorithm performance is illustrated in Table 5.40 and Figure 5.40. 143 Figure 5.36: The lookahead acceleration algorithm performance for the letter “G” movie with linear intensity increase input 144 Time Baseline QL Lookahead QL (TB PR =5) Lookahead QL (TB PR =10) 1 11807 11807 11807 2 6641 6641 6641 3 7559 7559 7559 4 7264 7242 7242 5 6338 6535 6535 6 5981 6260 6260 7 5719 5881 5881 8 5792 5498 5493 9 7017 6404 5368 10 9218 8519 5422 11 9257 9135 3568 12 3299 3410 2938 13 3149 3160 3025 14 3249 3259 3121 15 3357 3362 3217 Table 5.37: The lookahead acceleration algorithm queue lengths for the letter “G” movie with non-linear intensity increase input Time Bucket Speed Up Error TB PR =5 9.60% 17% TB PR =10 10.60% 45% Table 5.38: The lookahead acceleration algorithm performance for the letter “G” movie with non-linear intensity increase input 145 Figure 5.37: The lookahead acceleration algorithm queue lengths for the letter “G” movie with non-linear intensity increase input 146 Figure 5.38: The lookahead acceleration algorithm performance for the letter “G” movie with non-linear intensity increase input 147 Time Baseline QL Lookahead QL (TB PR =5) Lookahead QL (TB PR =10) 1 11746 11746 11746 2 6659 6659 6659 3 9278 9278 9278 4 9131 9106 9106 5 10539 10700 10700 6 11685 12068 12068 7 11697 11956 11956 8 11945 11690 11705 9 12408 11316 11593 10 12956 11477 11957 11 11442 10031 10078 12 3155 3268 3707 13 2992 2992 2992 14 3088 3088 3088 15 3184 3184 3184 Table 5.39: The lookahead acceleration algorithm queue lengths for the moving letter “G” movie with linear intensity increase input Time Bucket Speed Up Error TB PR =5 9.60% 23% TB PR =10 7.30% 30% Table 5.40: The lookahead acceleration algorithm performance for the moving letter “G” movie with linear intensity increase input 148 Figure 5.39: The lookahead acceleration algorithm queue lengths for the moving letter “G” movie with linear intensity increase input 149 Figure 5.40: The lookahead acceleration algorithm performance for the moving letter “G” movie with linear intensity increase input 150 Time Baseline QL SR QL (0.92Vmax) SR QL (0.90Vmax) SR QL (0.88Vmax) 1 2367 2367 2367 2367 2 4750 4750 4750 4750 3 5339 3327 4660 5216 4 2799 2207 2707 2806 5 2290 1723 2251 2292 6 2306 1958 2206 2305 7 2395 2154 2369 2394 Table 5.41: The self-regulating acceleration algorithm queue lengths for the random dots input Version Speed Up Error 0.92Vmax 16.90% 21% 0.90Vmax 4.21% 12% 0.88Vmax 0.53% 2% Table 5.42: The self-regulating acceleration algorithm performance for the random dots input 5.4 Advanced Inputs 5.4.1 Random Dots Self-Regulating Event queue growth is illustrated in Table 5.41 and Figure 5.41. Algorithm performance is illustrated in Table 5.42 and Figure 5.42. Event-Merging Event queue growth is illustrated in Table 5.43 and Figure 5.43. Algorithm performance is illustrated in Table 5.44 and Figure 5.44. 151 Figure 5.41: The self-regulating acceleration algorithm queue lengths for the random dots input 152 Figure 5.42: The self-regulating acceleration algorithm performance for the random dots input 153 Time Baseline QL Event-Merging QL (TB HC =1) Event-Merging QL (TB HC =2) 1 2367 2367 2367 2 4750 4750 4750 3 73 73 73 4 2270 2139 2139 5 1630 1316 1328 6 1601 1284 1495 7 1113 950 1011 8 1178 999 1024 9 1003 823 864 10 738 598 660 11 276 530 539 Table 5.43: The event-merging acceleration algorithm queue lengths for the random dots input Figure 5.43: The event-merging acceleration algorithm queue lengths for the random dots input 154 TB HC Outputs Erroneous Outputs Highly Erroneous % Mean Error Error Speed Up 1 317991 63266 78% 0.36Vmax 19% 37.37% 2 317991 73762 83% 0.40Vmax 23% 42.1% Table 5.44: The event-merging acceleration algorithm performance for the random dots input 155 Figure 5.44: The event-merging acceleration algorithm performance for the random dots input 5.4.2 Tommy Trojan Head We carried out simulations for the input illustrated in Figure 5.45. The baseline simula- tion outputs are illustrated in Figure 5.46. Event-Merging Event queue growth is illustrated in Table 5.45 and Figure 5.47. Algorithm performance is illustrated in Table 5.46 and Figure 5.48. 156 Figure 5.45: Tommy Trojan head Figure 5.46: Images represent Input, photoreceptor layer output, horizontal-cell layer output and bipolar-cell layer output Time Baseline QL Event-Merging QL (TB HC =2) Event-Merging QL (TB HC =3) 1 13404 13404 13404 2 23221 23221 23221 3 6297 6297 6297 4 18203 13300 13281 5 4510 3053 3056 6 3285 2509 2087 7 2456 2098 1913 8 1923 1777 1606 9 1569 1507 1419 10 980 850 862 11 467 1127 1069 12 241 829 820 13 173 1293 1193 14 61 1161 1138 15 0 1391 1369 Table 5.45: The event-merging acceleration algorithm queue lengths for the tommy trojan head input 157 Figure 5.47: The event-merging acceleration algorithm queue lengths for the tommy trojan head input 158 TB HC Outputs Erroneous Outputs Highly Erroneous % Mean Error Error Speed Up 2 4586570 926208 83% 0.37Vmax 20% 49.20% 3 4586570 956955 84% 0.39Vmax 20% 52.60% Table 5.46: The event-merging acceleration algorithm performance for the tommy trojan head input 159 Figure 5.48: The event-merging acceleration algorithm performance for the tommy trojan head input Self-Regulating Event queue growth is illustrated in Table 5.47 and Figure 5.49. Algorithm performance is illustrated in Table 5.48 and Figure 5.50. 160 Time Baseline QL SR QL (0.92Vmax) SR QL (0.90Vmax) SR QL (0.88Vmax) 1 13404 13404 13404 13404 2 23221 23221 23221 23221 3 41381 37811 39084 41381 4 4696 4427 4595 4694 5 4229 3843 4122 4208 Table 5.47: The self-regulating acceleration algorithm queue lengths for the tommy trojan head input Figure 5.49: The self-regulating acceleration algorithm queue lengths for the tommy trojan head input 161 HC STOP Outputs Erroneous Outputs Highly Erroneous % Mean Error Error Speed Up 0.88Vmax 663841 693 45% 0.13Vmax 0% 0% 0.90Vmax 663841 25064 58% 0.25Vmax 3% 2.9% 0.92Vmax 663841 34581 61% 0.20Vmax 5% 4.9% Table 5.48: The self-regulating acceleration algorithm performance for the tommy trojan head input 162 Figure 5.50: The self-regulating acceleration algorithm performance for the tommy trojan head input 163 5.5 Dynamic Time-Bucket Value In Section 6.2.1, we proposed the following formula to compute the dynamic bucket: TB i =k(aC iv + b C ir ) The equation basically implies that the time bucket inversely proportional to slope, and directly proportional to the integral. However, for this equation to work properly, the values of the constants must be chosen carefully and tuned inside the simulation. We propose another formula that can be executed in the simulation without the need to tune various constants: TB i = P i C ir . Cvmax Civ such that: P i is the smallest monotonic period of time in the neuronn i recent history. C vmax is the maximum produced volume by any neuron in the system. For the response represented in Figure 5.51,P i value is10. To compute the value for P i , the simulation continuously compares the new state with the previous state to detect monotonic periods and then determine the P i value. We chose the minimum period to reduce error, choosing higher values may increase the simulation speed but decrease the accuracy. Figure 5.51: A neuron response 164 The factor Cvmax Civ role is to assign higher time-bucket value for nodes that produce larger integral as they may influence the simulation more, and we are willing to spend more simulation time on them. 165 Chapter 6 Conclusions and Future Work 6.1 Conclusions 1 - The lookahead acceleration algorithm performs poorly for chaotic responses that do not follow any pattern since it is hard to carry out any form of extrapolation on these responses. 2 - The event-merging acceleration algorithm performs well for nodes receiving many events within short periods of time since this increases the probability of merging larger number of events. 3 - Self-regulating algorithm works well for slowly changing responses since it reduces the mean number of reported events from those nodes. 4 - Each acceleration algorithm has strengths and weaknesses given the various param- eters of the simulation experiment including the inputs and outputs of each node. Thus, it is important to dynamically select the appropriate algorithm for each node to increase the simulation acceleration algorithm performance. We discuss dynamic algorithm selection in Section 6.2.1. 5 - Simulation configuration decisions must be educated. To create a sound reti- nal model that produces the expected biological outputs, the model configuration including delays and connections must be biologically plausible. For example, in 166 Rubtina’s baseline model the delay between the photoreceptor and the bipolar cell is equal to the summation of the delay between the photoreceptor and the hori- zontal cell and the delay between the horizontal cell and the bipolar cell. These delay values will allow the processed outputs from both the photoreceptor and the horizontal cell to arrive at the bipolar cell at the same time. 6 - The time-bucket value has a big impact on the simulation performance. Choosing the right time bucket for each node is a key to increase the simulation performance. As the inputs and outputs of the simulation nodes vary during the simulation, the time bucket must be dynamically changed. 7 - The lookahead acceleration algorithm must utilize several extrapolation methods including piecewise linear and polynomial extrapolation. This will allow the algo- rithm to increase the simulation performance for larger spectrum of responses. The recent node responses can be used to select the appropriate extrapolation method. 8 - The acceleration simulation algorithm error relative to the baseline is measured by both the percentage of differences between the two simulations and the mean error value. Depending on the percentage of differences only did not give a clear idea for the self-regulating algorithm performance. The self-regulating algorithm ignores low-value events, resulting in larger number of differences between the accelerated simulation and the baseline. Upon closer inspection of the self-regulating sim- ulation algorithm output for various inputs, it appears that the output important characteristics like edges were not influenced as the algorithm ignore low-value events. Thus, considering the mean error value in the error analysis is important to understand the true algorithm performance. 167 9 - Rubtina is a cellular-level simulation since each cell is an autonomous system. The cellular-level implementation offers modelers higher flexibility to test various retinal and neural models. 6.2 Future Work 6.2.1 Dynamic Simulation Acceleration Initiatives Dynamic simulation acceleration attempts to modify the acceleration algorithms described in Section 3.4 during simulation to enable these algorithms to make more intel- ligent decisions in order to enhance their performance. Dynamic Time Bucket Computation The time-bucket value depends on each node’s response. We propose two factors to measure each node’s sensitivity to the inputs. Once response sensitivity is quantified, we can determine a suitable value for the time bucket. Some simulations have homogeneous nodes that will have the same time bucket across the entire simulation. However, when the simulation’s nodes are heterogeneous, the lookahead algorithm will compute time- bucket values for each node type. Given a node n i that receives a set of events E in the simulation future, we can measure the sensitivity of the node by measuring the response rate of change, and the produced activity volume of the node. The rate of change for noden i , C ir , can be approx- imated using the following equation: C ir = 1 M P M1 j=1 jf(E j )f(E j1 )j Such that f is the transfer function for the noden i , and M is the number of events that 168 will be considered from the noden i ’s response history. A node with slow response may have a bigger time bucket than a node that changes its response quickly. Figure 6.1 illustrates two responses with different rates of change such that the node with the response on the left could have a bigger time bucket than the node with the response to the right. Figure 6.1: The rate-of-change factor illustration. The second factor is the volume of the transfer function produced byn i , and this is approximated using the following equation: C iv = 1 M P M1 j=1 f(E j )+f(E j1 ) 2 Such that f is the transfer function for the noden i , and M is the number of events that will be considered from the noden i response history. The value ofM represents that number of samples that are considered from the node’s 169 history to predict its future. The modeler is expected to set an initial value forM. Figure 6.2 illustrates two responses with different response volume. Figure 6.2: The produced-volume factor illustration The time bucket will be inversely proportional to C ir and proportional to C iv values for each node. Thus, the time-bucket value will be:k(aC iv + b C ir ) =TB i Such that: a controls the importance of the node’s volume b controls the importance of the node’s rate of change k is a heuristic factor. The values for a, b and k are determined by the modeler. The time-bucket value for a certain node will be dynamically updated given the node’s recent history. Figure 6.3 illustrates the dynamic time bucket assignment. The time- bucket value for the neuron n i is updated whenever a set of M events has been processed by n i . The value of M can be initially set by the modeler and it can be modified during 170 simulation to more suitable value. The flowchart illustrated a variable R that is incre- mented whenever the neuron n i receives an event, and whenever it receives M events, it is ready to carry out an update of the time-bucket value. Figure 6.3: Dynamic time-bucket assignment Intelligent Self-Regulating Algorithm It would be useful to enhance the self-regulating algorithm to develop a more intelligent way to ignore low-value events. The self-regulating algorithm ignores an event when the percentage of change compared to the previous event is small or the event potential is below certain small threshold. However, the definition of a low-value event could be 171 advanced to include the influence of the event on the receiving neuron’s post-synaptic neurons. For example, if an event may make big change on its destination neuron, how- ever, this big change does not make any change on that neuron’s post-synaptic neurons, then we may consider this event a low-value event. Given a node n k sending events to a set of nodes S=n 1 , n 2 , .., n m as demonstrated in Figure 6.4. Each node in S can track it’s sensitivity to the input from node K. If most of the nodes in S are not sensitive to n k , then, the simulation will reduce processing node n k events as they are not influencing the simulation. Implementing this idea require adding a “presynaptic sensitivity table” to track the influence of the presynaptic neuron’s inputs on the neuron. Figure 6.4: Several connected neurons Event Merging Frequency The frequency of performing event-merging operations can be adjusted dynamically. The simulation may carry out these operations whenever a new event is inserted into the sec- ondary queue resulting in a large number of secondary queue scans. However, events are going to be merged when there is a set of events that falls within a time bucket, so it is not 172 Merging Frequency Event/Meta Event Ratio Computation Cost Speed Up Error Low Higher Lower Lower Higher High Lower Higher Higher Lower Table 6.1: Event merging trade-offs necessary to carry out frequent secondary queue scans before the secondary queue accu- mulates a sufficient number of events that can be merged together. The event-merging frequency trade-offs with algorithm computation cost and performance are illustrated in Table 6.1. We propose tracking the event mean inter-arrival time for each destination neuron for the recent period of simulation time. This period of time could be several time- bucket time units. The maximum number of events that can be merged into one meta event is the multiplication of the mean inter-arrival time by the length of the time bucket. As illustrated in Table 6.1, creating meta events with the highest number of events will produce the worst performance but the lowest computation cost. However, creating meta events with the smallest number of events (2 events) will produce the highest performance and highest computation cost. The number of events for each destination can be implemented using a lookup table with the destination address as a key, and the value stored is the number of events. This table must be update whenever we add or remove events from the secondary queue. Using the maximum number of events in a time bucket, and tracking the number of events in a secondary queue for each neuron, a decision can be made regarding the event-merging frequency by the modeler, given his simulation goals and available hardware. 173 Furthermore, the frequency of event-merging can be tuned dynamically given the simulation current activity. Such that, the simulation will attempt to increase the event- merging frequency during low-activity periods of the simulation (e.g. utilizing inactive simulation threads), and reduce the activity during high-activity periods. Dynamic Algorithm Selection Each acceleration algorithm performed differently given the inputs and outputs of the simulation node. The lookahead acceleration algorithm might be the best choice dur- ing a certain period of time for a particular node, but the self-regulating algorithm is better choice during another period of time for the same node. Thus, we need a way to change the simulation acceleration algorithm dynamically for every node to decide which algorithm will be enabled for a particular node and for which period of time. Thus, a “selector” algorithm is needed to enable and disable various acceleration algorithms dynamically for every simulation node. The selector algorithm may use the recent node responses and the consistency of these responses to select an acceleration algorithm. For example, if the recent responses for a particular node are slowly changing, the selector algorithm may choose the self-regulating algorithm. 6.2.2 Retinal Models In Rubtina, we have carried out our retinal models using simple models. Our models did not capture high biological details that other modelers may desire. However, Rubtina is implemented in a way to allow the modeler to add new models in simple way as explained in Section 4.5.2. In Section 6.2.2, we describe quite more detailed photoreceptor model that may be desired by some modelers, the described model can implemented and will run inside Rubtina. Other models can be formulated for other retinal cells. 174 Photoreceptor Cell Micro-Model In the absence of light, the photoreceptor response is stable at40mv. When the photore- ceptor is exposed to a light flash, the photoreceptor hyperpolarizes. That is, the membrane potential becomes more negative because of the reduction of positive ions inside the pho- toreceptor. The hyperpolarization leads to lower Ca2 + concentration, which results in the activation of a regulatory mechanism called light adaptation that regulates the membrane potential by bringing it back to the photoreceptor resting potential. Similar photoreceptor responses have been recorded by many researchers [5, 3, 4]. Figure 6.5: Biological retinal photoreceptor response [76] Photoreceptors absorb light photons through a complicated biochemical process con- verting the light into graded potentials. The human retina contains over 120 million rod cells and 6 million cone cells. Rods function mainly in dim light and provide black-and- white vision, while cones support daytime vision and the perception of color, making them more important for retinal prosthetic designers. The photoreceptors are able to absorb light since they contain visual pigments. When light falls into the photoreceptor, a certain fraction of the light photons will be absorbed by the photoreceptor depending on the cone outer-segment width, and the cone wavelength sensitivity. Photon absorp- tion will cause cone hyperpolarization (the membrane potential becomes more negative). 175 Once the photopigments are bleached and the cone receives the horizontal cell feedback, the cone membrane potential slowly will depolarize. Figure 6.6 illustrates the cone con- nectivity diagram. Figure 6.7 presents the cone-processing flowchart. Figure 6.6: Cone connectivity to other retinal cell in the human’s retina (C: cone and H: horizontal cell) The Michaelis-Menten saturation function is used to model the retinal cones. Figure 6.8 represents the Michaelis-Menten dynamics, which is represented by the equation: V = VmaxS km+S Such that: V : reaction rate 176 Figure 6.7: Cone processing flowchart V max : maximum reaction rate K m : the concentration S value at which the reaction rate is half ofV max S : substrate concentration 177 Figure 6.8: Michaelis-Menten saturation function [72] For a cone, the Michael-Menten saturation function is given using the following formula [98]: Vm(I) Vmax = I I+(kr+Ia kr k b ) Such that: V m : cone membrane potential V max : cone maximum membrane potential I: light intensity I half : light intensity that produces response withVmax=2 K r : half saturation factor K b : half pigment bleach constant I a : spatio-temporal ambient illumination intensity, modeled using the following formula: I a (x;y;t)=c 1 I B (x;y;t)+c 2 I S (x;y;t) Where: c 1 andc 2 : constants I B (x;y;t): bleached light intensity 178 I S (x;y;t): horizontal cells feedback intensity The c 1 and c 2 constants are used to determine the contribution of I B and I S . The bleach light intensity is determined by the photoreceptor bleaching time constant and can be computed by convolving the light intensity with a low-pass filter with a time constant equal to the photoreceptor bleaching time. Thus, bleaching light intensity is modeled using the following formula: I B (x;y;t)=I(x;y;t)K(t; bleach ) Such that: K(t; bleach ) is a low pass filter. bleach : bleaching time constant The Horizontal cells feedback intensity represents the feedback from horizontal cells to the photoreceptor that is modeled using the formula: I(x;y;t)G(x;y; horizontal )K(t; horizontal ) The parameter horizontal represents the width of the horizontal cell Gaussian filter, and horizontal is the horizontal cell time constant. Then, the cone membrane potential is: Vm(I) Vmax = I I+(kr+(c 1 I B (x;y;t)+c 2 I S (x;y;t)) kr k b ) And, the final cone model is: C(x;y;t)=V m (x;y;t)G(x;y; cone )K(t; cone )H(x;y;t) Such that: G(x;y; cone ): gaussian operator of width cone H(x;y;t): horizontal cell feedback 179 Bibliography [1] Analog VLSI Implementation of Neural Systems (The Springer International Series in Engineering and Computer Science). 1st ed. Springer, 1989. ISBN: 0792390407. URL: http : / / www . amazon . com / exec / obidos / redirect ? tag = citeulike07-20&path=ASIN/0792390407. [2] D´ avid B´ alya et al. “A CNN framework for modeling parallel processing in a mammalian retina: Research Articles”. In: Int. J. Circuit Theory Appl. 30.2-3 (2002), pp. 363–393. ISSN: 0098-9886. DOI: 10.1002/cta.v30:2/3. URL: http: //dx.doi.org/10.1002/cta.v30:2/3. [3] D. A. Baylor and M. G. F. Fuortes. “Electrical responses of single cones in the retina of the turtle”. In: The Journal of Physiology 207.1 (1970), pp. 77–92. ISSN: 0022-3751. URL: http://jp.physoc.org/content/207/1/77.abstract. [4] D. A. Baylor, M. G. F. Fuortes, and P. M. O’Bryan. “Receptive fields of cones in the retina of the turtle”. In: The Journal of Physiology 214.2 (1971), pp. 265–294. ISSN: 0022-3751. URL: http://jp.physoc.org/content/214/2/265.abstract. [5] D. A. Baylor, A. L. Hodgkin, and T. D. Lamb. “The electrical response of turtle cones to flashes and steps of light”. In: The Journal of Physiology 242.3 (1974), pp. 685–727. URL: http://jp.physoc.org/content/242/3/685.abstract. [6] David Beeman. Introduction to Realistic Neural Modeling. http://www.brains- minds-media.org/archive/218. 2009 (accessed November 10, 2009). [7] Computational Neurobiology Lab of The Salk Institute for Biology Study. MCell. http://www.mcell.cnl.salk.edu/. 2009 (accessed October 23, 2009). [8] Kim T. Blackwell. “Modeling Calcium Concentration and Biochemical Reac- tions”. In: Brains, Minds and Media 2 (2005). 180 [9] National Federation of the Blind. Blindness statistics. http://www.nfb.org/nfb/ blindness statistics.asp. 2009 (accessed August 26, 2009). [10] Kwabena A. Boahen and Andreas G. Andreou. “A contrast sensitive silicon retina with reciprocal synapses”. In: Advances in neural information processing sys- tems. V ol. 4. 1992, pp. 764–772. URL: http://citeseerx.ist.psu.edu/viewdoc/ summary?doi=10.1.1.43.8246. [11] James M. Bower and David Beeman. The book of GENESIS (2nd ed.): exploring realistic neural models with the GEneral NEural SImulation System. New York, NY , USA: Springer-Verlag New York, Inc., 1998. ISBN: 0-387-94938-0. URL: http://portal.acm.org/citation.cfm?id=275590. [12] J. K. Bowmaker and H. J. Dartnall. “Visual pigments of rods and cones in a human retina.” In: The Journal of physiology 298 (1980), pp. 501–511. ISSN: 0022-3751. URL: http://jp.physoc.org/cgi/content/abstract/298/1/501. [13] Dr. Thomas Caceci. The Retinal Tunic. http : / / education . vetmed . vt . edu / curriculum/VM8054/EYE/RETINA.HTM. 2009 (accessed August 6, 2009). [14] Srihari Cadambi, Chandra S. Mulpuri, and Pranav N. Ashar. “A Fast, Inexpensive and Scalable Hardware Acceleration Technique for Functional Simulation”. In: Design Automation Conference 0 (2002), pp. 570+. DOI: 10.1109/DAC.2002. 1012690. URL: http://dx.doi.org/10.1109/DAC.2002.1012690. [15] Cadence. Cadence Incisive Palladium series. http://www.cadence.com/products/ sd/palladium series/pages/default.aspx. 2010 (accessed August 2nd, 2010). [16] CalC The Calcium Calculator. Victor Matveev. http://www.calciumcalculator. org/. 2009 (accessed October 23, 2009). [17] Nicholas T. Carnevale and Michael L. Hines. The NEURON Book. Cambridge University Press, 2006. ISBN: 0521843219. URL: http://www.amazon.com/exec/ obidos/redirect?tag=citeulike07-20&path=ASIN/0521843219. [18] Shannon Carpentier, Maria Knaus, and Miyoung Suh. “Associations between Lutein, Zeaxanthin, and Age-Related Macular Degeneration: An Overview”. In: Critical Reviews in Food Science and Nutrition 49.4 (2009), pp. 313–326. ISSN: 1040-8398. DOI: 10.1080/10408390802066979. URL: http://dx.doi.org/10.1080/ 10408390802066979. 181 [19] R. Carrillo et al. “Event-driven simulation of neural population synchronization facilitated by electrical coupling”. In: Biosystems 87.2-3 (2007), pp. 275–280. ISSN: 03032647. DOI: 10.1016/j.biosystems.2006.09.023. URL: http://dx.doi.org/ 10.1016/j.biosystems.2006.09.023. [20] K. M. Chandy and J. Misra. “Asynchronous distributed simulation via a sequence of parallel computations”. In: Commun. ACM 24.4 (1981), pp. 198–206. ISSN: 0001-0782. DOI: 10.1145/358598.358613. URL: http://dx.doi.org/10.1145/ 358598.358613. [21] Sue-Yeon Choi et al. “Encoding Light Intensity by the Cone Photoreceptor Synapse”. In: Neuron 48.4 (2005), pp. 555–562. ISSN: 08966273. DOI: 10.1016/ j.neuron.2005.09.011. URL: http://dx.doi.org/10.1016/j.neuron.2005.09.011. [22] A. Y . Chow et al. “The artificial silicon retina microchip for the treatment of vision loss from retinitis pigmentosa.” In: Archives of ophthalmology 122.4 (2004), pp. 460–469. ISSN: 0003-9950. DOI: 10.1001/archopht.122.4.460. URL: http://dx.doi.org/10.1001/archopht.122.4.460. [23] C. A. Curcio et al. “Human receptor topography.” In: The Journal of comparative neurology 292.4 (1990), pp. 497–523. ISSN: 0021-9967. DOI: 10.1002/cne.9029 20402. URL: http://dx.doi.org/10.1002/cne.902920402. [24] A. Delorme et al. “SpikeNET: A simulator for modeling large networks of inte- grate and fire neurons”. In: Neurocomputing 26-27 (1999), pp. 989–996. ISSN: 09252312. DOI: 10.1016/S0925-2312(99)00095-8. URL: http://dx.doi.org/10. 1016/S0925-2312(99)00095-8. [25] Arnaud Delorme and Simon J. Thorpe. “SpikeNET: an event-driven simulation package for modelling large networks of spiking neurons.” In: Network (Bristol, England) 14.4 (2003), pp. 613–627. ISSN: 0954-898X. URL: http://view.ncbi. nlm.nih.gov/pubmed/14653495. [26] E Deschutter. “Nodus, a user-friendly neuron simulator for Macintosh comput- ers”. In: Neural Systems: Analysis and Modeling (1993), pp. 113–119. [27] RMagick developers. RMagick. http://rmagick.rubyforge.org/. 2011 (accessed August 22, 2011). 182 [28] Dawei W. Dong and Joseph J. Atick. “Temporal Decorrelation: A Theory of Lagged and Nonlagged Responses in the Lateral Geniculate Nucleus”. In: Net- work. V ol. 6. 1995, pp. 159–178. URL: http://citeseerx.ist.psu.edu/viewdoc/ summary?doi=10.1.1.47.808. [29] Dr and Joachim K. Anlauf. Distributed, Event Driven Simulation of Spiking Neu- ral Networks Dipl.-Ing. Cyprian Graßmann. URL: http://citeseerx.ist.psu.edu/ viewdoc/summary?doi=10.1.1.91.7939. [30] R. Eckmiller. “Learning retina implants with epiretinal contacts.” In: Ophthalmic research 29.5 (1997), pp. 281–289. ISSN: 0030-3747. URL: http://view.ncbi.nlm. nih.gov/pubmed/9323719. [31] C. Enroth-Cugell and J. G. Robson. “The contrast sensitivity of retinal ganglion cells of the cat.” In: The Journal of physiology 187.3 (1966), pp. 517–552. ISSN: 0022-3751. URL: http://view.ncbi.nlm.nih.gov/pubmed/16783910. [32] J. Evans. “Antioxidant supplements to prevent or slow down the progression of AMD: a systematic review and meta-analysis.” In: Eye (London, England) 22.6 (2008), pp. 751–760. ISSN: 0950-222X. DOI: 10.1038/eye.2008.100. URL: http: //dx.doi.org/10.1038/eye.2008.100. [33] J. R. Evans and K. Henshaw. “Antioxidant vitamin and mineral supplements for preventing age-related macular degeneration.” In: Cochrane database of sys- tematic reviews (Online) 1 (2008). ISSN: 1469-493X. DOI: 10.1002/14651858. CD000253.pub2. URL: http://dx.doi.org/10.1002/14651858.CD000253.pub2. [34] John H. Byrne Douglas A. Baxter Evyatar Av-Ron Michael J. Byrne. “SNNAP: A Tool for Teaching Neuroscience”. In: Brains, Minds & Media 1 (2008). [35] External Limiting Membrance. http://en.wikipedia.org/wiki/External limiting membrane. Wikipedia, 2009 (accessed August 8, 2009). [36] Eye. http://en.wikipedia.org/wiki/Eye. Wikipedia, 2009 (accessed August 16, 2009). [37] R. G. Foster et al. “Circadian photoreception in the retinally degenerate mouse (rd/rd)”. In: Journal of Comparative Physiology A: Neuroethology, Sensory, Neural, and Behavioral Physiology 169.1 (1991), pp. 39–50. DOI: 10 . 1007 / BF00198171. URL: http://dx.doi.org/10.1007/BF00198171. 183 [38] R. M. Fujimoto. “Parallel discrete event simulation”. In: WSC ’89: Proceedings of the 21st conference on Winter simulation. Washington, D.C., United States: ACM, 1989, pp. 19–28. ISBN: 0-911801-58-8. DOI: 10.1145/76738.76741. URL: http://dx.doi.org/10.1145/76738.76741. [39] Joel S. Glaser and Alfredo A. Sadun. Neuro-Ophthalmology. Third Edition. Lip- pincott Williams & Wilkins, 1999. Chap. 4, pp. 75–92. ISBN: 0781717299. URL: http://www.amazon.com/exec/obidos/redirect?tag=citeulike07-20&path=ASIN/ 0781717299. [40] M. A. Goodale and A. D. Milner. “Separate visual pathways for perception and action.” In: Trends in neurosciences 15.1 (1992), pp. 20–25. ISSN: 0166-2236. URL: http://view.ncbi.nlm.nih.gov/pubmed/1374953. [41] Dan Goodman and Romain Brette. “Brian: a simulator for spiking neural net- works in Python”. In: BMC Neuroscience 9.Suppl 1 (2008), P92+. ISSN: 1471- 2202. DOI: 10.1186/1471-2202-9-S1-P92. URL: http://dx.doi.org/10.1186/1471- 2202-9-S1-P92. [42] Grace. Grace. http://plasma-gate.weizmann.ac.il/Grace/. 2009 (accessed October 23, 2009). [43] NEURON Group. NEURON. http : / / www . neuron . yale . edu / neuron/. 2009 (accessed November 14, 2009). [44] XuDong Guan and Hui Wei. “Realistic Simulation on Retina Photoreceptor Layer”. In: Hainan Island, China, Apr. 2009, pp. 179–184. DOI: 10.1109/JCAI. 2009.148. URL: http://dx.doi.org/10.1109/JCAI.2009.148. [45] Lance Hahn. Compartmental Models. http://retina.anatomy.upenn.edu/ lance/ modelmath/model compartment.html. 2009 (accessed October 23, 2009). [46] Deborah E. Hannula, Daniel J. Simons, and Neal J. Cohen. “Imaging implicit perception: promise and pitfalls”. In: Nat Rev Neurosci 6.3 (2005), pp. 247–255. ISSN: 1471-003X. DOI: 10.1038/nrn1630. URL: http://dx.doi.org/10.1038/nrn16 30. [47] Dyonne T. Hartong, Eliot L. Berson, and Thaddeus P. Dryja. “Retinitis pigmen- tosa.” In: Lancet 368.9549 (2006), pp. 1795–1809. ISSN: 1474-547X. DOI: 10. 1016/S0140-6736(06)69740-7. URL: http://dx.doi.org/10.1016/S0140-6736(06) 69740-7. 184 [48] D. Helbing. “Micro- and macro-simulation of freeway traffic”. In: Mathematical and Computer Modelling 35.5-6 (2002), pp. 517–547. ISSN: 08957177. DOI: 10. 1016/S0895-7177(02)80019-X. URL: http://dx.doi.org/10.1016/S0895-7177(02) 80019-X. [49] Richard F. Spaide Howard F. Fine. Exudative Age-Related Macular Degenera- tion: Causes and Treatment. http://cme.medscape.com/viewarticle/545207 3. 2009 (accessed August 27, 2009). [50] David H. Hubel. Eye, Brain and Vision. http://hubel.med.harvard.edu/book/b9. htm. 2009 (accessed August 23, 2009). [51] David H. Hubel. Eye, Brain and Vision. http://hubel.med.harvard.edu/book/b8. htm. 2009 (accessed August 24, 2009). [52] David H. Hubel. Eye, Brain and Vision. http://hubel.med.harvard.edu/book/b12. htm. 2009 (accessed August 27, 2009). [53] David H. Hubel. Eye, Brain and Vision. http://hubel.med.harvard.edu/book/b13. htm. 2009 (accessed August 27, 2009). [54] jgarridoalcazar. Edlut. http://code.google.com/p/edlut/. 2010 (accessed March 29, 2010). [55] J. B. Jonas, U. Schneider, and G. O. Naumann. “Count and density of human retinal photoreceptors.” In: Graefe’s archive for clinical and experimental oph- thalmology = Albrecht von Graefes Archiv f¨ ur klinische und experimentelle Oph- thalmologie 230.6 (1992), pp. 505–510. ISSN: 0721-832X. URL: http://view.ncbi. nlm.nih.gov/pubmed/1427131. [56] P. A. Juliano and E. Rosenbaum. “A novel SCR macromodel for ESD circuit simulation”. In: Washington, DC, USA, pp. 14.3.1–14.3.4. DOI: 10.1109/IEDM. 2001.979499. URL: http://dx.doi.org/10.1109/IEDM.2001.979499. [57] Arie Kaufman. “A micro-macro simulation model for a signalized network”. In: ANSS ’80: Proceedings of the 13th annual symposium on Simulation. Piscataway, NJ, USA: IEEE Press, 1980, pp. 87–101. URL: http://portal.acm.org/citation.cfm? id=808869. 185 [58] Robert K. Koenekoop et al. “Novel RPGR mutations with distinct retinitis pig- mentosa phenotypes in French-Canadian families.” In: American journal of oph- thalmology 136.4 (2003), pp. 678–687. ISSN: 0002-9394. URL: http://view.ncbi. nlm.nih.gov/pubmed/14516808. [59] Y . C. Liang et al. “A neural-network-based method of model reduction for the dynamic simulation of MEMS”. In: Journal of Micromechanics and Microengi- neering 11.3 (2001), pp. 226+. ISSN: 0960-1317. DOI: 10.1088/0960-1317/11/3/ 311. URL: http://dx.doi.org/10.1088/0960-1317/11/3/311. [60] Collin J. Lobb et al. “Parallel Event-Driven Neural Network Simulations Using the Hodgkin-Huxley Neuron Model”. In: PADS ’05: Proceedings of the 19th Workshop on Principles of Advanced and Distributed Simulation. Washington, DC, USA: IEEE Computer Society, 2005, pp. 16–25. ISBN: 0-7695-2383-8. DOI: 10.1109/PADS.2005.18. URL: http://dx.doi.org/10.1109/PADS.2005.18. [61] Textensor Ltd. PSICS - the Parallel Stochastic Ion Channel Simulator. http:// www.psics.org/index.html. 2009 (accessed October 23, 2009). [62] B. D. Lubachevsky. “Efficient distributed event-driven simulations of multiple- loop networks”. In: Commun. ACM 32.1 (1989), pp. 111–123. ISSN: 0001-0782. DOI: 10.1145/63238.63247. URL: http://dx.doi.org/10.1145/63238.63247. [63] Macular degeneration. http : / / en . wikipedia . org / wiki / Macular degeneration. Wikipedia, 2009 (accessed August 17, 2009). [64] M. A. Mahowald. “Silicon retina with adaptive receptors”. In: Society of Photo- Optical Instrumentation Engineers (SPIE) Conference Series. Ed. by B. P. Mathur and C. Koch. V ol. 1473. Society of Photo-Optical Instrumentation Engi- neers (SPIE) Conference Series. 1991, pp. 52–58. URL: http://adsabs.harvard. edu/cgi-bin/nph-bib query?bibcode=1991SPIE.1473...52M. [65] Misha Mahowald. “VLSI analogs of neuronal visual processing : a synthesis of form and function”. PhD thesis. 1992. URL: http://www.worldcat.org/isbn/. [66] E. Margalit et al. “Retinal Prosthesis for the Blind”. In: Survey of Ophthalmology 47.4 (2002), pp. 335–356. ISSN: 00396257. DOI: 10.1016/S0039-6257(02)0031 1-9. URL: http://dx.doi.org/10.1016/S0039-6257(02)00311-9. 186 [67] D. Marr and T. Poggio. “A Computational Theory of Human Stereo Vision”. In: Proceedings of the Royal Society of London. Series B, Biological Sciences 204.1156 (1979), pp. 301–328. ISSN: 00804649. DOI: 10.2307/77509. URL: http: //dx.doi.org/10.2307/77509. [68] Victor Matveev, Arthur Sherman, and Robert S. Zucker. “New and corrected sim- ulations of synaptic facilitation.” In: Biophysical journal 83.3 (2002), pp. 1368– 1373. ISSN: 0006-3495. DOI: 10.1016/S0006-3495(02)73907-6. URL: http://dx. doi.org/10.1016/S0006-3495(02)73907-6. [69] P. A. McNaughton. “Light response of vertebrate photoreceptors.” In: Physio- logical reviews 70.3 (July 1990), pp. 847–883. ISSN: 0031-9333. URL: http:// physrev.physiology.org/content/70/3/847.abstract. [70] Carver Mead. Analog VLSI and Neural Systems. 1st. Addison Wesley Publishing Company, 1989, pp. 257–278. ISBN: 0201059924. URL: http://www.amazon. com/exec/obidos/redirect?tag=citeulike07-20&path=ASIN/0201059924. [71] Inc. Merck Co. Age-Related Macular Degeneration: Retinal Disorders: Merck Manual Home Edition. http://www.merck.com/mmhe/sec20/ch234/ch234b.html. 2009 (accessed August 8, 2009). [72] MichaelisMenten kinetics. http : / / en . wikipedia . org / wiki / Michaelis - Menten kinetics. Wikipedia, 2011 (accessed August 10, 2011). [73] Jayadev Misra. “Distributed discrete-event simulation”. In: ACM Comput. Surv. 18.1 (1986), pp. 39–65. ISSN: 0360-0300. DOI: 10.1145/6462.6485. URL: http: //dx.doi.org/10.1145/6462.6485. [74] Laney T. Mobley. Comparison of Graded Potentials and Action Potentials. http: //kcfac.kilgore.cc.tx.us/mobleypageap2/unit1/actionandgradedpotentials.htm. 2009 (accessed August 25, 2009). [75] NRCAM. Virtual Cell. http://www.nrcam.uchc.edu/. 2009 (accessed October 23, 2009). [76] OpenWetWare. Phototransduction. http : / / openwetware . org / wiki / BIO254 : Phototransduction. 2010 (accessed April 16 2010). [77] G. Osterberg. “Topography of the layer of rods and cones in the human retina”. In: Acta Ophthalmol. 13 (1935), pp. 1–102. 187 [78] Alice C. Parker and Adi N. Azar. “A hierarchical artificial retina architecture”. In: ed. by ´ Angel B. Rodr’iguez V´ azquez. V ol. 7365. 1. Dresden, Germany: SPIE, 2009, pp. 736503+. URL: http://scitation.aip.org/getabs/servlet/GetabsServlet? prog=normal&id=PSISDG007365000001736503000001&idtype=cvips&gifs= yes. [79] Anjali Patel. A General Overview on Visual Perception. http : / / serendip . brynmawr.edu/bb/neuro/neuro00/web1/Patel.html. 2009 (accessed Dec 25, 2009). [80] Jim Perlewitz. Computational Neuroscience Software. http://home.earthlink.net/ perlewitz/sftwr.html. 2009 (accessed October 23, 2009). [81] “Petroleum Reservoir Simulations: A Basic, Approach, J. H. Abou-Kassem, S. M. Farouq Ali, and M. R. Islam, Gulf Publishing Company, Houston, Texas, USA, 2006, 480 pages, ISBN: 0-9765113-6-3”. In: Petroleum Science and Tech- nology 25.6 (2007), pp. 825–826. DOI: 10.1080/10916460701435008. URL: http: //dx.doi.org/10.1080/10916460701435008. [82] M. Pi´ orkowski et al. “TraNS: realistic joint traffic and network simulator for V ANETs”. In: SIGMOBILE Mob. Comput. Commun. Rev. 12.1 (2008), pp. 31– 33. ISSN: 1559-1662. DOI: 10.1145/1374512.1374522. URL: http://dx.doi.org/ 10.1145/1374512.1374522. [83] G. F. Poggio and T. Poggio. “The Analysis of Stereopsis”. In: Annual Review of Neuroscience 7.1 (1984), pp. 379–412. DOI: 10.1146/annurev.ne.07.030184. 002115. URL: http://dx.doi.org/10.1146/annurev.ne.07.030184.002115. [84] Receptive field. http://en.wikipedia.org/wiki/Receptive field. Wikipedia, 2010 (accessed March 31, 2010). [85] Pittsburgh Supercomputing Center Researchers. PGENESIS. http://www.psc. edu/Packages/PGENESIS/. 2010 (accessed March 31, 2010). [86] Retinal Degenerations: Biology, Diagnostics, and Therapeutics (Ophthalmology Research). 1st ed. Humana Press, 2007. ISBN: 1588296202. URL: http://www. amazon.com/exec/obidos/redirect?tag=citeulike07-20&path=ASIN/1588296202 . [87] Retinal Ganglion Cell. http : / / en . wikipedia . org / wiki / Retinal ganglion cell. Wikipedia, 2009 (accessed August 18, 2009). 188 [88] Retinitis pigmentosa. http : / / en . wikipedia . org / wiki / Retinitis pigmentosa. Wikipedia, 2009 (accessed August 17, 2009). [89] R. W. Rodieck. “Quantitative analysis of cat retinal ganglion cell response to visual stimuli.” In: Vision research 5.11 (1965), pp. 583–601. ISSN: 0042-6989. URL: http://view.ncbi.nlm.nih.gov/pubmed/5862581. [90] Austin Roorda and David R. Williams. “The arrangement of the three cone classes in the living human eye”. In: Nature 397.6719 (Feb. 1999), pp. 520–522. ISSN: 0028-0836. DOI: 10.1038/17383. URL: http://dx.doi.org/10.1038/17383. [91] Eduardo Ros et al. “Event-driven simulation scheme for spiking neural networks using lookup tables to characterize neuronal dynamics”. In: Neural Comput. 18.12 (2006), pp. 2959–2993. ISSN: 0899-7667. DOI: 10.1162/neco.2006.18. 12.2959. URL: http://dx.doi.org/10.1162/neco.2006.18.12.2959. [92] Michael H. Rowe. “Trichromatic Color Vision in Primates”. In: News in Phys- iological Sciences 17.3 (2002), pp. 93–98. URL: http : / / physiologyonline . physiology.org/cgi/content/full/17/3/93. [93] W. Rushton and G. Henry. “Bleaching and regeneration of cone pigments in man”. In: Vision Research 8.6 (June 1968), pp. 617–631. ISSN: 00426989. DOI: 10.1016/0042-6989(68)90040-0. URL: http://dx.doi.org/10.1016/0042-6989(68) 90040-0. [94] Scholarpedia. Retina. http://www.scholarpedia.org/article/Retina. 2009 (accessed August 25, 2009). [95] E. De Schutter. A consumer guide to neuronal modeling software. http://www. tnb.ua.ac.be/publications/pub002/TNB pub2.shtml. 2009 (accessed October 26, 2009). [96] ScienceDaily. New Parallelization Technique Boosts Computers’ Ability to Model Biological Systems. http : / / www. sciencedaily. com / releases / 2011 / 06 / 110609112911.htm. 2011 (accessed November 26, 2011). [97] Samir Shah et al. Information Processing in Primate Retinal Cone Pathways: A Model. 1993. URL: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.32. 9070. 189 [98] Samir Shah et al. Information Processing in Primate Retinal Cone Pathways: Experiments and Results. 1993. URL: http : / / citeseerx . ist . psu . edu / viewdoc / summary?doi=10.1.1.32.6909. [99] R. G. Smith. “NeuronC: a computational language for investigating functional architecture of neural circuits.” In: Journal of neuroscience methods 43.2-3 (1992), pp. 83–108. ISSN: 0165-0270. URL: http : / / view. ncbi . nlm . nih . gov / pubmed/1405746. [100] Robert G. Smith. Retsim. http://retina.anatomy.upenn.edu/ rob/neuronc.html. 2009 (accessed November 12, 2009). [101] Laboratory of Neural Engineering of the University of Southern California. EONS. http://www.synaptic-modeling.com. 2009 (accessed October 23, 2009). [102] Peter Sterling and Jonathan B. Demb. The Synaptic Organization of the Brain. 5th ed. Oxford University Press, USA, 2003. Chap. 6, pp. 217–269. ISBN: 019515956X. URL: http : / / www . amazon . com / exec / obidos / redirect ? tag = citeulike07-20&path=ASIN/019515956X. [103] R. Newton A. Sangiovanni-Vincentelli T. Quarles D. Pederson and Christopher Wayne. SPICE. http : / / bwrc . eecs . berkeley. edu / classes / icbook / spice/. 2010 (accessed April 27th, 2010). [104] Blue Brain Team. Frequently Asked Questions. http://bluebrain.epfl.ch/page1892 4.html. 2010 (accessed August 3rd, 2010). [105] Cadence Team. Accelerating Analog Simulation With Full Spice Accuracy. http: //www.cadence.com/rl/Resources/white papers/spectre turbo wp.pdf. 2010 (accessed August 6th, 2010). [106] Cadence Team. Virtuoso Spectre Circuit Simulator. http://www.cadence.com/ products/rf/spectre circuit/pages/default.aspx. 2010 (accessed August 6th, 2010). [107] DOE Artificial Retina Project Team. About the Artificial Retina Project. http: //www.artificialretina.energy.gov/about.shtml. 2009 (accessed July 21, 2009). [108] Eve team. Simulation Acceleration With EVE’s ZeBu. http://www.eve-team.com/ solutions/simulation-acceleration.php. 2010 (accessed August 2nd, 2010). 190 [109] NS2 Team. The Calendar Queue Scheduler. http://www.isi.edu/nsnam/ns/doc/ node35.html. 2010 (accessed August 3rd, 2010). [110] NS2 Team. The Network Simulator - ns-2. http://www.isi.edu/nsnam/ns/. 2010 (accessed August 6th, 2010). [111] Odysse research team. Virtual Retina - Customization. http://www-sop.inria.fr/ odyssee/software/virtualretina/customization.shtml. 2009 (accessed September 17, 2009). [112] PSICS Team. PSICS Examples. http : / / www. psics . org / examples . html. 2009 (accessed November 9, 2009). [113] verisity. Xtreme Server. http://www.verisity.com/products/xtreme server.html. 2010 (accessed August 2nd, 2010). [114] Visual phototransduction. http : / / en . wikipedia . org / wiki / Visual phototransduction. Wikipedia, 2009 (accessed Dec 25, 2009). [115] C. Visweswariah and R. A. Rohrer. “Piecewise approximate circuit simulation”. In: IEEE Transactions on Computer-Aided Design of Integrated Circuits and Sys- tems 10.7 (1991), pp. 861–870. ISSN: 02780070. DOI: 10.1109/43.87597. URL: http://dx.doi.org/10.1109/43.87597. [116] James D. Weiland, Wentai Liu, and Mark S. Humayun. “RETINAL PROSTHE- SIS”. In: Annual Review of Biomedical Engineering 7.1 (2005), pp. 361–401. DOI: 10.1146/annurev.bioeng.7.060804.100435. URL: http://dx.doi.org/10.1146/ annurev.bioeng.7.060804.100435. [117] Adrien Wohrer and Pierre Kornprobst. “Virtual Retina : A biological retina model and simulator, with contrast gain control”. In: Journal of Computational Neuro- science 26.2 (2009), pp. 219–249. DOI: 10.1007/s10827-008-0108-4. URL: http: //dx.doi.org/10.1007/s10827-008-0108-4. [118] Kwoon Y . Wong, Felice A. Dunn, and David M. Berson. “Photoreceptor Adap- tation in Intrinsically Photosensitive Retinal Ganglion Cells”. In: Neuron 48.6 (2005), pp. 1001–1010. ISSN: 08966273. DOI: 10.1016/j.neuron.2005.11.016. URL: http://dx.doi.org/10.1016/j.neuron.2005.11.016. 191 [119] J. L. Wyatt, D. L. Standley, and W. Yang. “The MIT Vision Chip Project: ana- log VLSI systems for fast image acquisition and early vision processing”. In: Robotics and Automation, 1991. Proceedings., 1991 IEEE International Confer- ence on. 1991, 1330–1335 vol.2. DOI: 10.1109/ROBOT.1991.131797. URL: http: //dx.doi.org/10.1109/ROBOT.1991.131797. [120] Kareem Zaghloul. “A Silicon Implementation of a Novel Model for Retinal Pro- cessing”. PhD thesis. University of Pennsylvania, 2001. 192 Appendix A: The Anatomy of the Mammalian Retina: an Overview The Visual Pathway Vision is a complicated process that requires the mammalian eye and brain to work together [40]. The visual pathway in most mammalians goes through the eye, optic nerve, optic chiasma, optic tract, lateral geniculate body (LGN), optic radiation, and visual cor- tex. The eye is often compared to a camera. In that analogy, the retina would be analogous to the camera film, the pupil would be analogous to the aperture, and the cornea and lens would be analogous to the camera lens. Light falling onto the eye is refracted 1 as it passes through the cornea. The cornea is the front part of the mammalian’s eye that covers the iris, pupil, and the anterior chamber. It plays the initial role in light refraction, which ensures that visual scene image is focused regardless of the object distance. It then passes through the pupil 2 and is further refracted by the lens. The diameter of the pupil 1 Refraction is the wave directional change. 2 The pupil is a circular black opening located in the center of the eye’s iris that allows light to reach the retina. 193 is controlled by the iris 3 . The cornea and lens act together as a compound lens to project an inverted focused image onto the retina. In the human retina, over 125 million retinal photoreceptors [77] transduce light into action potentials, and the remaining retinal cells (horizontal, bipolar, amacrine and ganglion cells) are responsible for “processing” these action potentials to deliver information efficiently to the brain. Figure 9 illustrates the anatomy of the eye. Figure 9: The anatomy of the human eye [36] The retinal output is sent to the brain through the optic nerve. The optic nerve contains the axons of the retinal ganglion cells (about 1 million axons in the human retina). The optic nerves from both eyes meet and cross at the optic chiasm. The visual information for both eyes is combined and then splits based on the visual field (right or left). Both the right and left visual fields are sent to the left and right halves of the mammalian brain 3 The iris is a membrane that is responsible for controlling the diameter of the pupil, which determines the amount of light reaching the retina. 194 respectively. As such, the right side of primary visual cortex is responsible for the left half of the visual scene, and the left side of the primary visual cortex is responsible for the right half of the visual scene. For visual scene information to reach the corresponding visual primary cortex, it passes through the optic tract. Both left, and right optic tracts terminate in the lateral geniculate nucleus (LGN). LGN has six layers that are mapped to different ganglion cell types. The exact function of the LGN is not well known. The authors in [28] believe that LGN accomplishes temporal decorrelation. The neurons of the LGN relay the visual scene information to the primary visual cortex (V1). V1 is the front end of the visual cortex, it is located at the back of the brain. Then, the information flows through remaining visual cortex areas including V2, V3, V4 and area V5. The visual cortex functions include motion detection, object location, and face recognition [40]. Figure 10 represents the visual pathway. The Mammalian Retina Overview The retina is a light sensitive neural tissue partially covering the inner surface of the eye. It is considered a part of the brain itself ( 0.5% of the brain’s weight) [102]. The main role of the retina is to transduce light into action potentials that carry the visual scene information in an efficient way. The retina consists of several complex neural layers that contains the five major retinal cells including: photoreceptors (rods and cones), horizontal cells, bipolar cells, amacrine cells, and ganglion cells. The ganglion cells lie in the innermost part of the retina, the photoreceptors cells, and horizontal cells lie in the outermost layer and the bipolar and amacrine cells lie in the middle. The photoreceptor cells perform light transduction, while other carry out various visual scene processing on 195 Figure 10: The visual pathway [46] the transduced signals. The human retina has over 125 million photoreceptors, and this helps create quality vision. Those 125 million photoreceptors eventually create input to 1.2 million ganglion cells. Thus, each ganglion cell indirectly receives inputs from 100 photoreceptors, which is an indication on how much the retinal cells encode the visual scene information to deliver an efficient visual image to the brain. The visual scene information path in the retina is classified into direct and indirect paths [50]. The direct path is the path from photoreceptors to bipolar cells to ganglion cells. It is direct since the information from these cells flows forward toward the ganglion cells. The indirect path includes horizontal cells feeding information back to photorecep- tors and/or bipolar cells, and amacrine cells feeding information to ganglion cells and/or 196 bipolar cells. Both horizontal and amacrine cells perform lateral inhibition on the other cells. The photoreceptors, horizontal and bipolar cells (distal neurons) respond to stimulus with a sustained graded potential, while ganglion and some amacrine cells (proximal neurons) respond using action potentials (all or none). There are two possibilities for the distal neurons to use graded potentials [94]. They do not need to carry information long distances like other cells, and the graded potential give these cells the ability to give more information than the action potential. The action potential amplitude usually is 0v or 100mv, while the graded potential varies from 1mv to 50mv [74]. When light falls onto the retina, it passes through the entire thickness of the retina before reaching the photoreceptors because of the retinal cells arrangement. The pho- toreceptors transduce light into graded potential and the remaining retinal cells carry out various visual scene processing. There is not complete understanding of what processing is carried out in the retina. However, the visual processing includes luminance adaptation, contrast enhancement, edge detection and possibly more. The retinal output is carried on the ganglion cells axons (optic nerve), and will be sent to the primary visual cortex. The retinal structure is quite backward in the sense that photoreceptors are located on the very back of retina rather than the front end. Hubel [51] gave two possible jus- tifications for this curious backward structure. The first is the existence of Melanin in the back of the retina, which absorbs the light preventing it from being reflected back that may create unclear vision. The second possible reason is the photoreceptors need to recreate the light sensitive material and this process requires the Melanin, which is in the outermost part of the retina. One of the most pigmented areas in the retina is the macula. The center of the macula is the fovea, which is the most sensitive area in the retina to light and responsible for 197 our sharp central vision. It is very sensitive to light since it is densely packed with cones (199,000 cones/mm 2 [23]). The central retina extends around the fovea for about 6 mm and then the peripheral retina. The edge of the retina is defined by the ora serrata. Another area in the retina is the optic disk. It is completely free of photoreceptors, which results in the lack of vision in the associated area in the visual field. The optic disk is often called “the blind spot”. The optic disk is the area where ganglion cell axons exit the eye to form the optic nerve. It appears as a vertical oval white area with average dimensions of 1.76 mm horizontally by 1.92 mm vertically [39]. Retinal Layers The retina contains ten different layers (illustrated in Figure 11) including [13] Figure 11: The ten retinal layers [13] 198 1 - The Pigment Epithelium Layer (PE) is a single layer of cuboidal cells. Its main visual task is to increase the contrast of the visual scene image since it absorbs the light that would otherwise be reflected back inwards towards the retinal pho- toreceptors. Furthermore, this layer stores Vitamin A that is produced by R&CL layer. 2 - The Rods and Cones Layer (R&CL) contains part of the photoreceptors; the actual sensitive part of the photoreceptor cell. The sensitive part consists mainly of rhodopsin, which is a pigment that immediately photobleaches when exposed to light. There is a continuous rhodopsin cycling processing (formation/breakdown) that happens due to the photobleaching process of these pigments. 3 - The Outer Limiting Membrane (OLM) isolates the inner layers of the retina from any potential harmful material in the blood circulation, forming a blood-retina bar- rier. Figure 12 illustrates the OLM layer position just near the base of a photore- ceptor. The OLM separates the photosensitive part of the photoreceptors from the remaining part, which ensures that all materials formed during the cyclic formation and breakdown of rhodopsin will not reach other parts in the retina. 4 - The Outer Nuclear Layer (ONL) contains the nuclei and cell bodies of the photore- ceptor cells. 5 - The Outer Plexiform Layer (OPL) contains numerous synapses between the pho- toreceptors and integrator neurons. These integrating elements include horizontal, and bipolar cells. 6 - The Inner Nuclear Layer (INL) contains cell bodies and nuclei of the bipolar and horizontal cells. 199 Figure 12: The outer limiting membrane layer position [35] 7 - The Inner Plexiform Layer (IPL) contains the synapses of both bipolar and gan- glion cells. 8 - The Ganglion Cell Layer (GCL) contains the cell bodies of the ganglion cells. 9 - The Nerve Fiber Layer (NFL) contains the axons of the ganglion cells that will send the action potentials to the visual cortex through the optic nerve. 10 - The Inner Limiting “Membrane” (ILM) seals the retina from the potential harmful materials in the vitreous chamber (the space between the lens and the retina). Retinal Disorders Millions of people suffer from various eye disorders. The National Federation of the Blind (NFB) [9] reports that 1.3 Million people are legally blind in the US. NFB antic- ipates that this number going to be around 1.6 million in 2015. There are quite large 200 number of eye disorders [9] that affect various parts in the eye including sclera, cornea, iris, lens, choroid and retina. Retinal disorders include retinal detachment (the retina detaches from the choroid), Retinoschisis (retinal layers separate from each other), Mac- ular edema (swollen macula cause problem in central vision), retinitis pigmentosa (RP), and age-related macular degeneration (AMD). Age-Related Macular Degeneration (AMD) AMD causes progressive damage to one of the most vital areas in the retina; the macula, resulting in gradual irreversible loss of the central vision [71]. There are two types of AMD, wet and dry. The dry form is more common, however, the wet form causes the majority of blindness among AMD patients. In the dry form, cellular debris accumulates between the retina and the choroid, and the retina becomes eventually detached. Both eyes may be affected in the dry form at the same time. In the wet form, blood vessels (hence the description “wet”) grow from the choroid behind the retina, and the retina becomes eventually detached just like the dry form. The wet form develops in one eye first, but eventually may affect both eyes. Some of the AMD causes include the following: 1 - High blood pressure 2 - Older age 3 - Smoking 4 - Fair skin: Patients with very fair skin has demonstrated increased risk [49]. How- ever, in [86] the authors question that fair skin may increase AMD risks. 5 - Genetic abnormalities 6 - A diet low in omega-3 fatty acids. 201 Wet macular degeneration starts off as the dry type. In wet macular degeneration, abnormal blood vessels develop under the macula. These vessels may leak fluid and blood under the retina. Eventually, a scar tissue develops under the retina. AMD treatment is quite limited. However, Some evidence supports a reduction in the risk of AMD with increasing the patient’s daily intake of lutein and zeaxanthin [18]. Other research studies [33, 32] could not show any relationship between antioxidants, and minerals and AMD. Retinitis Pigmentosa (RP) Retinitis Pigmentosa is a genetic retinal disorder that causes abnormalities of the photoreceptors or the retinal pigment epithelium layer [47]. The symptoms of RP start by night blindness followed by the loss of peripheral vision (no loss in central vision). Most RP patients do not become legally blind until their 40s or 50s [58]. RP patients see black surrounding in their visual field around the center of the image as illustrated in the left image in Figure 13. The RP treatment is quite limited, however, there are new promising products from ScyFIX MCN. ScyFIX [107] develops eye therapy products using micro-current neuro- modulation (MCN). ScyFIX MCN is the first approved device therapy in the world for treating Retinitis Pigmentosa (RP). Figure 13 illustrates the way AMD and RP patients perceive an image. Retinal Cells Photoreceptors A photoreceptor is a retinal neuron that is capable of phototransduction. The photore- ceptors absorb light photons and through a complicated biochemical process, convert the light into graded potential. The human retina contains about 120 million rod cells and 202 Figure 13: The figure illustrates an image (left) and the way it is perceived by AMD (center) and RP (right) patients. Images are taken from [88, 63] and merged together. 6 million cone cells. For hundreds of years, scientists thought that there were only two types of photoreceptor cells: cones and rods. However, recent research results found a third type in the ganglion cells [37, 118]. However, those cells are not involved directly in image forming function like the rods and cones. Rods function mainly in dim light and provide black-and-white vision, while cones support daytime vision and the perception of color. In the human’s retina, there are three major types of cones: red, green, and blue cones. Those three colors have different wavelengths, and the correspondent cone can only absorb that specific color wavelength. When light contains a mixture of red and green colors, only the red and green cones create a graded potential. Indeed, all colors that humans are able to perceive are a mixture of those three colors. People who suffer from color blindness have a problem with one or more cone types that cause a deficiency in perceiving colors. Rods function well in low light vision (during night). Figure 14 represents the absorption of rods, and various cones to light wavelengths. Photoreceptors respond to the wavelength and intensity of the visual scene. To deter- mine color, the visual system compares responses across the population of the red, green and blue cones [92]. The photoreceptors are able to absorb light since they contain visual 203 Figure 14: Photoreceptors absorption rates. Curve labeled as follows: 498 the absorbance of rods, 420 absorbance of blue cones, 534 absorbance of green cones, and 564 is the absorbance of the red cones [12]. pigments (rhodopsin). Furthermore, photoreceptors respond to various stimuli intensity with different graded potentials [21]. Horizontal Cells Horizontal cells are the laterally-interconnecting neurons in the OPL. Scientists believe that they inhibit the photoreceptors, or bipolar cells or both [52]. Once photoreceptors hyperpolarize, they hyperpolarize the horizontal cells that will depolarize the photore- ceptors. Hyperpolarization is a change in the potential of a cell’s membrane that makes it more negative. Once the photoreceptor hyperpolarizes, it reduces the release of its excitatory neurotransmitter (glutamate), then horizontal cells reduce the release of their inhibitory neurotransmitter (GABA) that has an inhibitory effect on the photoreceptors, 204 which leads to photoreceptor depolarization. This feedback from horizontal cells to pho- toreceptors appears to be a way to control the brightness, which may indicate that hori- zontal cells have a role in luminance adaptation. Horizontal cells possibly connect to the surround of the bipolar cells. Horizontal cells communicate among each other using gap junctions. A gap junction is an intercellular connection between cells that allows various molecules and ions to pass freely between cells. This communication allows certain horizontal cell outputs to be influenced by distal horizontal cells. They communicate with each other to create a spatially averaged map of light intensities. Bipolar Cells Bipolar cells are located between the photoreceptors and ganglion cells (hence the description, “bipolar”), they transmit the processed photoreceptor output to the ganglion cells. Each bipolar cell receives input from either cone cells or rod cells. As such, bipolar cells are designated as rod bipolar or cone bipolar cells based on their input source (rod or cone respectively). There are about ten forms of cone bipolar cells, and only one rod bipolar cell. The bipolar cells are center-surround structures in which the cell has a central area, and surround area. The output of the cell is a function of the inputs on the center and surround. The bipolar cell’s center mainly receives inputs from retinal photoreceptors. A possible input for the surround is the horizontal cells output. The cone bipolar cells can be classified into ON and OFF cells, based on their reactions to the photoreceptors’ inputs. The ON cell will depolarize when the photoreceptor hyperpolarizes, and vice versa for the OFF cell. The center-surround behavior appears to be a way to create edge detection as explained in Section 6.2.2. 205 Amacrine Cells Amacrine cells operate at the inner plexiform layer (IPL). There are a large number of amacrine cells types, most lacking axons. Amacrine cells work laterally affecting the bipolar cells by inhibiting them. Each amacrine cell type connects with a particular bipolar cell type, and generally has a particular type of neurotransmitter. These cells link the bipolar and ganglion cells together, creating an alternative indirect path between those cells. Most amacrine cells functions are probably unknown [53]. It is fair to say that out knowledge about these cells is quite limited. Ganglion Cells The ganglion cells are the back end of the retina. Their axons deliver the output to the visual cortex. They are center-surround structures just like the bipolar cells. The response profile of these cells is well described by a difference of Gaussian filter that is the basis for many edge detection algorithms used by computer scientists. Ganglion cells are also classified based on their spatial summation type. Cells showing linear spatial summation are termed X cells, and those showing non-linear summation are termed Y cells. The ganglion cells can be also classified based on their visual function [87]: 1 - M cells (aka Parasol): large center-surround receptive fields that are sensitive to depth, insensitive to color, and can respond to low-contrast stimuli. 2 - P cells (aka Midget): small center-surround receptive fields that are sensitive to color. 80% of ganglion cells are P cells. 3 - K cells (aka Bistratified): large center-only receptive fields that are sensitive to color, and moderate-contrast stimuli. 206 4 - Photosensitive ganglion cells [37, 118]: another photoreceptor cell that was discov- ered recently. Their functions are very different from the cones and rods as being non-image-forming photoreceptors. Their functions include controlling the pupil size, and synchronizing circadian rhythms 4 . 5 - Ganglion cells that controls the eye movement. Phototransduction in the Dark and Light The retinal photoreceptors are depolarized in dark, and hyperpolarized under bright light. Most neurons respond in a binary way, all-or-none (spike or do not spike) output. How- ever, this is not the case for retinal photoreceptors as their response is graded potential. Indeed all retinal cells respond using graded potentials except the amacrine and ganglion cells. In this section, we briefly discuss their operation in both dark and light conditions [114, 79]. In the Dark In the dark, high levels of cGMP 5 keep cGMP-gated sodium channels open, allowing a steady inward current, called the dark current. The dark current maintains neuron depolarization at -40 mV . The depolarization of the cell membrane opens voltage-gated calcium channels. An increased intracellular concentration of Ca2 + causes vesicles to release glutamate into the synaptic cleft. Glutamate release leads to the hyperpolariza- tion of the ON bipolar cells, and depolarization of the OFF bipolar cells. 4 Circadian rhythm is the daily biological cycle of living entities. 5 cGMP (Cyclic Guanosine Monophosphate) is a common regulator of ion channels. 207 Action/Cell Type ON Cell OFF Cell Light on center only Cell is firing rapidly Cell is not firing Light on surround only Cell is not firing Cell is firing rapidly No light at all Cell is not firing Cell is not firing Light on center and surround Cell is firing slowly Cell is firing slowly Table 2: The response of both the retinal ON and OFF ganglion cells for various stimuli focused on the center and/or surround of the retinal ganglion cells. In the Light When the light falls onto the photoreceptors, the photopigment molecules in the photore- ceptor outer segment absorb the photons. Due to the photon-absorbed energy, a series of chemical processes occurs on the photopigment molecules, which leads to cGMP con- centration reduction and therefore the cGMP-gated sodium channels will close. Clos- ing these channels will hyperpolarize the neuron, and the photoreceptor will reduce the amount of the released glutamate. Reducing the neurotransmitter will depolarize the ON bipolar cells and hyperpolarize the OFF bipolar cells. Retinal Spatial Encoding The retina spatially encodes the visual scene image to reduce the amount of information sent through the optic nerve. The retina has 100 times more Photoreceptor cells than ganglion cells that implies that the capacity of the optic nerve is limited and this is why encoding is important. These operations are carried out by the ON and OFF bipolar and ganglion cells; these cells often termed the center-surround structures. Table 2 represents the responses of the center-surround structures (both ON and OFF ganglion cells) for various stimuli scenarios. 208 The center-surround structure operation is equivalent to the edge detection algorithms that are used to extract the edges in a digital image. Figure 15 represents the ON gan- glion cells response (right), and OFF ganglion cell response (left). The figure illustrates Hubel’s description of the Stephen Kuffler retinal ganglion cell’s response recording experiment [51]. He describes the ganglion ON cells response when a stimulus is over its center as follows. If we were to hear the ganglion cell response using a loud speaker, we would hear an initial spontaneous firing and once light goes on, we would hear barrage of signals that sounded like a machine gun firing. If we look at Figure 15, the second response on the left part of the figure illustrates what Hubel said. When stimulus was applied to the center only, immediate intense pulses showed up. When stimulus stayed on, the cells kept firing at a lower rate. Figure 15: ON and OFF ganglion cells [50] responding to various center and surround stimuli. Summary Large numbers of people through out the world suffer from various retinal disorders that may lead to peripheral or central vision loss or irreversible blindness due to retinal disorder like AMD and RP. The lack of medical cure for common retinal disorders is 209 a big motivation to look closer into the retinal physiology. The various eye parts work together to focus light on the retina. The retina transduces and “processes” the light and send the output to the brain. The retina ignores the less important details in the visual scene as it tries to deliver “interesting” information like edges. 210 Appendix B: Retinal Prosthesis Models Even though our work focus is to create a retinal simulation platform, looking at reti- nal prosthesis modeling is inspiring, and it is useful to understand prosthetic designers’ needs. In this section, we briefly review few retinal prosthesis projects and discuss their retinal models. In [116, 78], retinal prosthesis projects have been surveyed. We can classify all the work that has been done in the field of artificial vision into two major approaches: extraocular and intraocular. Extraocular approaches usually involve systems that do not have any implants inside the eye. Cortical prosthesis devices are an example of extraocular approaches. Intraocular approaches include systems fully or partially implanted inside the eye. Intraocular approaches differ in being either epiretinal or subretinal [66]. In epiretinal designs, the prosthesis device will have parts inside and outside the eye. However, in subretinal designs, the device will be fully implanted inside the eye. VLSI Analogs of Neuronal Visual Processing Mahowald et al. proposed a visual system that consists of a silicon retina, stereo- correspondence chip and silicon optic nerve [65]. The silicon retina is responsible for light transduction and signal processing. The stereo-correspondence chip determines the 211 location of the object in depth using stereo-vision algorithms. The silicon optic nerve performs digital communication using a technique called address-event representation. The silicon retina consists of transducers and silicon neurons placed on a 2D array. The silicon retina reduces the bandwidth of information by subtracting the average inten- sity levels from the visual scene light intensity. The proposed silicon retina has a spatial averaging circuit that is a resistive network, in which each node is connected to six other nodes. This resistive network mimics the behavior of the electrical network formed by the gap junctions among neighboring horizontal cells. Just like a human retina, each node potential is going to be affected by the neighboring nodes since resistors connect each node to others. The silicon retina includes an inhibition feedback from the resistive network to the photoreceptors. This performs a luminance adaptation phenomena similar to the human retina. The silicon retina is fabricated on a single separate chip. Two retinal models were proposed: a feed forward [70] and a feed backward [64] silicon retina. The difference between these two silicon retinas is the existence of feedback from the hori- zontal network to the photoreceptor in the feed backward retina. The models attempt to model the OPL only. There is not any amacrine or ganglion cell functionalities included in these models. Figure 16 shows that silicon retina representation that consists of a resistive network and a single pixel element (inside the circular window). These pixels are tiled in a hexag- onal array. The output of the photoreceptor and the resistive network node potential are differentiated and the output is connected back to the resistive network. Both the poten- tial of the resistive network and the output of photoreceptors are differentiated, and their output represents the bipolar cell output 212 Figure 16: A silicon retina implementation [65] The silicon optic nerve is the chip responsible of transferring information from the silicon retina to the stereo-correspondence chip. The chip uses an asynchronous digi- tal multiplexing technique called address-event representation. In this communication protocol, the address of the source neuron is being sent on the wire to inform receiver neurons that this particular neuron is firing. This encoding scheme reducesN wires to (1 +log 2 N). The stereo-correspondence chip’s main functionality is to create a sen- sation of stereopsis in the visual system [83]. The chip receives two images from two silicon retinas through the optic nerve, and creates an image with a sensation of depth (stereopsis). Mahowald et al. proposed stereo-correspondence chip used a stereo-vision algorithm [1] that is based on a modified version of a previous stereo-matching algorithm [67]. 213 Artificial Silicon Retina (ASR) Artificial Silicon Retina is an example of a subretinal intraocular device [22]. Drs. Alan and Vincent Chow of the Optobionics company developed the ASR chip. It is 2 mm diameter and 25 m thick silicon-based subretinal microchip that has 5000 electrodes. Each electrode represents a pixel of the visual scene image. The chip uses the incident light to power itself. This passive device is composed of tiny solar cells that, when activated by light, will send an electrical signal through the secondary neurons of the retina down the optic nerve to the brain. A potential problem with such passive devices is that they may generate too little power to be effective. This problem is very prominent in dark areas as there is not enough light to power the chip. In [22] the ASR team reported that the chip was implanted in 6 patients. During a 6 to 18 month follow up period with all the patients; the ASR chips functioned electrically well. Visual improvements were reported, and no signs of implant rejection or infection showed up. The team suggested having larger clinical trial to further evaluate their chip safety and efficiency. U.S. Department of Energy (DOE) Artificial Retina Project One of the most successful and inspiring retinal prosthesis projects is the U.S. Depart- ment of Energy (DOE) Artificial Retina Project led by Dr. Mark Humayun of the Doheny Eye Institute at USC. It is a multi-institutional effort to develop an implantable micro- electronic epiretinal device that restores vision to people blinded by retinal diseases. The DOE team is working on three models of this chip. These models primarily differ in the number of electrical electrodes that the model has. The goal is to have a large number of electrodes since this will increase the visual scene resolution. 214 Model 1 of the chip has a 16-electrode array that allows the implanted electronic chip to wirelessly communicate with a camera mounted on a pair of glasses. Model 1 of their artificial retina was implanted in 6 patients whose ages ranged from 56 to 77 at the time of implant between 2002 and 2004. This implant now enables patients to detect when lights are on or off, describe an objects motion (right, left, up, and down), and count individual items. Patients have very limited vision with model 1 because the implant has only 16 electrodes. Model 2 has 60 electrodes, and it is underway in clinical trials; having more electrodes will help patients have better vision. Model 3 is currently under development. It will have about 200 electrodes. In [116], the DOE team stated that having at least 1000 electrodes (using simulations) in the implanted chip would allow blind patients to be able to perform quite advanced visual functions like facial recognition and reading. The DOE artificial retina consists of extraocular and intraocular systems. The extraocular systems consist of a small camera that is built into a pair of glasses, battery unit, and a visual processing unit. The intraoc- ular unit contains the array of electrodes. The extraocular camera captures the visual scene, and transmits it to the visual processing unit that processes the image. The image is transmitted to the intraocular unit using a telemetry system. The intraocular system is a stimulator that stimulates the remaining functional part of the retina. 215
Abstract (if available)
Abstract
The biological retina is a complex system with intricate neural processing, dense connectivity and massive number of neurons. Unlike most neurons, most retinal neurons respond with graded potentials that will increase the retinal system neural-state simulation update rate. The retinal neurons connect densely to each other using various connectivity patterns including feedforward, feedback and laterally leading to an even higher neural-state update rate. Considering also the massive number of neurons, modeling the retinal system is very challenging. ❧ The lack of medical cures for common retinal disorders has motivated many research groups to study the retina and to attempt to model this complicated system. Some of the retinal projects have created retinal prosthesis devices fabricated on chips, while other groups have created retinal simulators. Retinal simulators are used as powerful learning tools to understand the retinal neurons and their role in vision. For retinal prosthesis research groups, retinal simulators are important tools to determine a suitable quantity of neurons in the prosthetic device design (e.g. number of photoreceptors to allow face recognition). ❧ We describe an efficient event-driven retinal simulator that speeds up the simulation by detecting periods of ""inconsequential activity"" in the neural responses, and carries out approximations or eliminations of simulation events during these periods. We implemented several simulation acceleration algorithms including a lookahead event-driven algorithm that carries out the approximations and self-regulating and event-merging event-driven algorithms that eliminate events. ❧ We conducted various simulation experiments using a variety of inputs and showed the acceleration and error trade-offs for each input. The performance of each algorithm has been computed relative to the baseline event-driven simulation that does not include any acceleration algorithms. ❧ The lookahead event-driven algorithm carries out a piecewise linear approximation to approximate the future response of the node being simulated given its last two responses. In the lookahead event-driven algorithm we propose, events that are scheduled at a node Ni with a delivery time less than or equal to Tcurrent+TBi can be delivered to Ni at Tcurrent, such that TBi is Ni time bucket. The lookahead carries out an approximation of the time bucket events using the last two computed responses. Reducing computations will increase the simulation speed, but will possibly reduce accuracy. The time bucket value is determined by the modeler. ❧ Event merging and self-regulating aim to reduce the number of events in the system to increase the simulation speed while maintaining accurate results. Self-regulating operates to determine if the simulation can drop a low-value event without processing it. A low-value event is an event that does not influence the destination node's response or the influence is weak. The event-merging acceleration algorithm attempts to create more-complex meta events, each of which entails a set of events that occurs within a time bucket. The algorithm scans the event queue at certain times and attempts to create one event object (meta event) from several events that are scheduled to be delivered to the same destination with delivery times that fall in the same time bucket. ❧ The lookahead algorithm performed well for experiments with linearly changing responses. The self-regulating algorithm performed well for slowly changing responses since it reduces the mean number of reported events from the slowly changing nodes. And, the event-merging acceleration algorithm performed well for nodes receiving many events within short periods of time since this increases the probability of merging a larger number of events. ❧ Our main contributions are the following: First, we describe a comprehensive outer retinal simulator framework that is flexible and biomimetic, and processes realistic input, producing tangible output. Second, the simulator employs novel event-driven simulation acceleration strategies that are capable of increasing the simulation speed at much higher rate than the error increases. These strategies would be helpful on any continuous system simulation and are not unique to retinal or biological simulations.
Conceptually similar
PDF
Intraocular and extraocular cameras for retinal prostheses: effects of foveation by means of visual prosthesis simulation
PDF
Manipulation of RGCs response using different stimulation strategies for retinal prosthesis
PDF
Contributions to structural and functional retinal imaging via Fourier domain optical coherence tomography
PDF
Neuromorphic motion sensing circuits in a silicon retina
PDF
Towards a high resolution retinal implant
PDF
Improving stimulation strategies for epiretinal prostheses
PDF
Proactive detection of higher-order software design conflicts
PDF
A reference architecture for integrated self‐adaptive software environments
PDF
Dendritic computation and plasticity in neuromorphic circuits
PDF
Model-driven situational awareness in large-scale, complex systems
PDF
Ice: an impedance spectroscopy and atomistic simulation study
PDF
Simulating electrical stimulation and recording in a multi-scale model of the hippocampus
PDF
Therapeutic electrical stimulation strategies for neuroregeneration and neuroprotection of retinal neurons
PDF
Data-driven 3D hair digitization
PDF
Real-time simulation of hand anatomy using medical imaging
PDF
High-throughput methods for simulation and deep reinforcement learning
PDF
Sensing with sound: acoustic tomography and underwater sensor networks
PDF
A complex event processing framework for fast data management
PDF
Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells
PDF
A signal processing approach to robust jet engine fault detection and diagnosis
Asset Metadata
Creator
Azar, Adi N.
(author)
Core Title
Adaptive event-driven simulation strategies for accurate and high performance retinal simulation
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Computer Engineering
Publication Date
04/17/2012
Defense Date
01/13/2012
Publisher
Los Angeles, California
(original),
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
event driven,OAI-PMH Harvest,retina,retinal prosthesis,retinal simulation,simulation
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Parker, Alice C. (
committee chair
), Jenkins, Brian Keith (
committee member
), Medvidović, Nenad (
committee member
)
Creator Email
adi.azar@yahoo.com,adiazar@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-7121
Unique identifier
UC1110344
Identifier
usctheses-c3-7121 (legacy record id)
Legacy Identifier
etd-AzarAdiN-613.pdf
Dmrecord
7121
Document Type
Dissertation
Rights
Azar, Adi N.
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
event driven
retina
retinal prosthesis
retinal simulation
simulation