University of Southern California Dissertations and Theses

Adaptive octree meshing of the extracellular space in neural simulations

(USC Thesis Other)

Adaptive octree meshing of the extracellular space in neural simulations

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content ADAPTIVE OCTREE MESHING OF THE EXTRACELLULAR SPACE IN NEURAL SIMULATIONS by Christopher Benedict Charles Girard A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (BIOMEDICAL ENGINEERING) August 2023 Copyright 2023 Christopher Benedict Charles Girard Table of Contents List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Intracellular space: compartmental models . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Extracellular space: extending compartments . . . . . . . . . . . . . . . . . . . . 3 1.3.1 Analytical solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3.2 Implicit meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Extracellular space: mesh-based methods . . . . . . . . . . . . . . . . . . . . . . 5 1.4.1 Mesh generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4.2 Conductivity discretization . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4.3 Constraints & solving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.5 Potential improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Chapter 2: Equivalent-Circuit Transforms . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.2 Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 2d square: electrical equivalence is impossible . . . . . . . . . . . . . . . . 15 2.3.2 2d square: error minimization . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.3 Nodal interpolation without splitting . . . . . . . . . . . . . . . . . . . . . 18 2.4 Discussion and alternative formulations . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.1 Admittance method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.2 Mesh dual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.4.3 Trilinear finite element . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Chapter 3: Adaptive Octree Meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1 Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2.2 Octree-based meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 ii 3.2.3 Splitting metric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2.4 Netlist creation, boundary conditions, and solving . . . . . . . . . . . . . . 36 3.2.5 Validation tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.5.1 Common features . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.5.2 Conductance discretizations . . . . . . . . . . . . . . . . . . . . 39 3.2.5.3 Error metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.1 Discretization formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.2 Octree and uniform meshes . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4.1 Singularity effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4.2 Nonconforming elements . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4.3 Simplifications and shortcomings . . . . . . . . . . . . . . . . . . . . . . . 53 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Chapter 4: Electrode Representation and Activation Thresholds . . . . . . . . . . . . . 55 4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.2 Initial study: toy neuron . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.2.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.3 Realistic neuron model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.4 Discussion and conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Chapter 5: Dynamic Mesh Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.1 Single compartment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.1.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.1.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2 Toy multi-neuron model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.2.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3 Realistic model geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.4 Discussion and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Chapter 6: Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 A Size field parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 B Derivation of finite element method . . . . . . . . . . . . . . . . . . . . . . . . . 94 iii List of Figures 1.1 Examples of neurological devices in clinical and research applications . . . . . . . 2 1.2 Schematic of a basic neural compartment . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Schematic of an expanded compartmental model with explicit extracellular repre- sentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Overview of typical mesh-based simulation . . . . . . . . . . . . . . . . . . . . . 6 1.5 Simulated neural activity in a hippocampal slice following stimulation . . . . . . . 9 2.1 Characteristic impedances of a homogenous 2d square . . . . . . . . . . . . . . . . 13 2.2 Refinement of 2d square . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Solutions on 2d square . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4 Refinement error in split mesh by resistor values . . . . . . . . . . . . . . . . . . . 17 2.5 Refinement error at candidate minima of split mesh . . . . . . . . . . . . . . . . . 19 2.6 Schematic of diagonal refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.7 Refinement error in crossed mesh by resistor values . . . . . . . . . . . . . . . . . 21 2.8 Circuits resulting from different methods of discretizing conductivity in a cuboid region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.9 Schematic of admittance method voxel . . . . . . . . . . . . . . . . . . . . . . . . 24 2.10 Equivalent circuits of the mesh dual . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.11 Equivalent circuit of a trilinear finite element . . . . . . . . . . . . . . . . . . . . . 29 3.1 Illustration of meshing algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2 Diagram of octree structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.3 Solutions from admittance, FEM, and mesh dual methods . . . . . . . . . . . . . . 43 3.4 Accuracy of admittance, trilinear FEM, and dual methods applied to the same oc- tree mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 iv 3.5 Total computation time as a function of mesh size . . . . . . . . . . . . . . . . . . 45 3.6 Solution and error in XY plane (Z=0) for octree depth 10 . . . . . . . . . . . . . . 46 3.7 Error distribution as mesh resolution increases . . . . . . . . . . . . . . . . . . . . 47 3.8 Detail of error at non-source nodes in an octree mesh, maximum depth of 14 . . . . 48 3.9 Net error of uniform and octree meshes . . . . . . . . . . . . . . . . . . . . . . . . 49 3.10 Time requirements of simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.11 Overshoot as localℓ 0 approaches source radius . . . . . . . . . . . . . . . . . . . 51 4.1 Examples of extracellular potential at peak of stimulus . . . . . . . . . . . . . . . 58 4.2 Variation in activation threshold as a function of local mesh resolution . . . . . . . 59 4.3 Ratio of activation thresholds for meshed disk electrode vs. analytic point-source . 60 4.4 Extracellular potential of each compartment as mesh resolution varies . . . . . . . 61 4.5 Comparison of strength-duration curves for common models of myelinated neu- rons, reproduced from [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.6 Detail of MRG neuron model 25µm (left) and 150µm (right) away from the electrode. 63 4.7 Activation threshold of MRG neurons by mesh resolution and distance from electrode 64 5.1 Images of single-source simulation with dynamically adapted mesh . . . . . . . . . 69 5.2 Accuracy and computational speed of fixed and dynamically adapted meshes. . . . 70 5.3 Total error and computational time for fixed and dynamic meshes . . . . . . . . . . 71 5.4 Images of simulation with dynamically adapted mesh . . . . . . . . . . . . . . . . 74 5.5 Accuracy and computational speed of fixed and dynamically adapted meshes. . . . 75 5.6 Total error and computational time for fixed and dynamic meshes . . . . . . . . . . 76 5.7 Schematic of DBS electrode in hippocampus (a) and slice of minimal resolution mesh (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.8 Potential field in the hippocampus under pseudorandom stimulation. . . . . . . . . 80 5.9 Summary comparison of fixed and dynamically-adapted meshes in the hippocam- pal simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 v Abstract While compartmental models have scaled to millions of neurons across thousands of cores, high resource consumption and limited parallelism has restricted simulation of the extracellular domain to comparatively low-fidelity approximations. Moreover, neural activity exhibits significant spa- tiotemporal variation, so reallocating detail at each timestep could speed simulation and improve accuracy. This adaptation is considered as both an equivalent-circuit transform on the electrical network, and as an independent manipulation of mesh geometry. The former involves mutually incompatible constraints and is an unsuitable theoretical basis; the latter approach instead yields reasonable approximations. The result is a simple, easily parallelized, means of discretizing the extracellular space accord- ing to activity. By recursively splitting until a target size is reached, the entire mesh maintains an octree structure amenable to localized alteration. With element size directly proportional to dis- tance from the source, octree meshes achieve accuracy similar to a uniform grid of their smallest element, at a fraction of the computational burden. Point sources reveal a potential pitfall as element sizes approach the source’s radius, where potential in central elements significantly overshoots before additional mesh nodes enter the source and begin explicitly representing the boundary. Though naive simulation of injected current could thus yield inaccurate activation of nearby neural compartments, in practice the highly localized phenomenon does not affect overall excitability. vi Finally, dynamic mesh adaptation yields accuracy comparable to a single, fixed mesh over the course of local field potential and multi-electrode stimulus simulations, but is an order of magnitude faster. vii Chapter 1 Introduction 1.1 Motivations Aside perhaps for patch-clamp studies, which allow isolated access to a neuron’s surface or inte- rior, any investigation of a neuron’s electrical activity will be mediated by the extracellular space between it and the electrode. As figure 1.1 illustrates, these electrode-tissue systems arise in clin- ical and investigational use, in driving their activity and passively recording it, all in a variety of different anatomical regions. While nothing can truly replaceinvivo studies, they are difficult and costly to conduct—to say nothing of the added ethical concerns of animal experimentation. The ability to accurately simulate prior to confirming in reality is thus a major facilitator of research and development, especially as computational power continues to increase and information sharing becomes easier. Though technological progress does enable unprecedented size and scale of simulations, not all problems scale with the increase in performance. These neural-machine interfaces pose partic- ular challenges to modellers, as the two represent at very different geometric scales and dynamic behaviors, even when their interaction is governed by the same quantities of voltage and current. 1 (a) Hippocampal prosthetic (b) Retinal prosthetic (c) Peripheral nerve electrode array Figure 1.1: Examples of neurological devices in clinical and research applications, reproduced from [2, 3, 4], respectively. 1.2 Intracellular space: compartmental models Though their time-domain behavior is complex, the electrophysiology of neurons can be modeled well by elementary circuit theory. 19th century investigation into long-distance electrical trans- mission provided an analogy for neurons in cable theory; in this case, with transmission along the internal conductance influenced by the capacitance and resistance of the membrane. These contin- uous passive properties could then be approximated by discrete circuit components in a repeating pattern, effectively subdividing the length of the neuron (or transatlantic cable, in the original ap- plication). Following Hodgkin and Huxley’s characterization of ion channels, the electrical representation of the a membrane segment expanded to the circuit of figure 1.2[5]. Now, dynamic, nonlinear conductances describe the net behavior of the ion channels responsible for the action potential, in addition to the extant capacitance of cable theory. Though those particular dynamics are unique to neuron electrophysiology, the conceptual framework is still very much a product of circuit theory. 2 Figure 1.2: Schematic of a basic neural compartment 1.3 Extracellular space: extending compartments For studies focused on membrane dynamics alone, the extracellular space can be an afterthought; the neuron’s environment is highly conductive, so simply treating it as electrical ground is a conve- nient simplification. Of course, such an approach is useless if extracellular electrodes are in play; instead, the extracellular space around a compartment can be modeled by a similar arrangement of circuit elements, as illustrated in figure 1.3. Extracellular potential can now vary among compart- ments, possibly inducing an action potential; however, these potentials now need to be explicitly calculated to constrain the system. 1.3.1 Analytical solutions Most commonly, the extracellular space is modeled separately from the intracellular compartmen- tal models, with the results of one model applied as constraints to the other—extracellular voltage of the compartmental model, and transmembrane current as an extracellular source. Particularly when details of the neurons’ biophysics are the focus, it is common to simply use the analytical solutions of idealized cases rather than explicitly modelling the extracellular environment. 3 Figure 1.3: Schematic of an expanded compartmental model with explicit extracellular represen- tation Assuming a spherically-symmetric source 1 embedded in an infinite, homogenous, isotropically conductive medium, the potential field resulting from current injection is v(r)= i src 4πσr (1.1) where r is the Euclidean distance from the source’s center. Of course, no real scenario will actu- ally satisfy those assumptions, but at sufficient distance from inhomogeneities it is a reasonable approximation. This description has modeled both stimulating electrodes and spiking neural com- partments; alternatively, the latter can use the slightly more realistic line-source representation [1, 6, 7]. For a small number of simulated sources, it is relatively simple to determine the net potential field by the sum of their resulting fields, according to superposition. However, if there are many sources, or the field must be calculated at many points, even these relatively simple calculations can become burdensome. In a realistic scenario spanning multiple tissue types (and thus different electrical properties), these calculations may be misleading as well as inefficient. 1 including point sources 4 1.3.2 Implicit meshes Of course, this electrical network can be further extended by explicitly modelling conductances between extracellular nodes, linking compartments that are geometrically nearby but not topologi- cally connected. In instances of tight, fairly regular packing of cells, it can be feasible to calculate additional resistance terms between extracellular compartments of adjacent cells and simulate as a unified circuit. The bundled, parallel neurons of a nerve fiber allow easy calculation of the appro- priate conductance between neighbors, such as by Delaunay triangulation[8]. A similar approach has been employed in a smooth muscle model, where the uniformity of cell geometries result in a uniform 3d grid [9]. In practice, though, it is uncommon to simulate intracellular and extracellular space in this sort of unified model. A realistic arrangement of neurons will not necessarily provide exploitable geometric regularity, complicating the initial calculation of extracellular conductances. Moreover, their inclusion as explicit components of the network greatly increases the computational require- ments of the simulation [10, 11, 12, 13]. 1.4 Extracellular space: mesh-based methods Consequently, more detailed models typically rely on a mesh-based simulation such as the finite element method. A general description of these methods is illustrated in figure 1.4. 1.4.1 Mesh generation Fundamentally, the simulation domain is divided into many simple geometric regions, typically tetrahedra or hexahedra, known as elements; collectively, their arrangement is a mesh. Within each element, electrical properties are constant, and a simple mathematical expression (typically a low-order polynomial) describes the potential field. Thus, accuracy of the solution is generally tied to the resolution of the mesh: smaller elements better approximate the true field. 5 Discetize conductivity - Finite Element? - Admittance? - Graph Dual? Constrain boundaries - Merge equipotential nodes - Inject current at electrodes, compartments Mesh geometry - Conforming? - Structured? - Hex/tet/mixed? S O L U T I O N S E T U P Figure 1.4: A typical mesh-based simulation. A) Bulk electrical properties of tissues and elec- trodes are defined, as well as their geometric boundaries. B) The domain geometry is divided into a collection (the mesh) of geometric primitives (the elements); the type of mesh is categorized by the shapes employed, as well as the connectivity’s regularity (conforming or non-conforming) and representation in memory (structured or unstructured). Continuous conductivity across each ele- ment is approximated by discrete conductances/resistances, as a function of the element geometry and regional conductivity. D) Equipotential nodes are combined, and the fixed potentials and cur- rents applied to the electrical network. E) Nodal analysis yields a sparse system of equations that is solved for all nodal potentials. Values of the field between nodes can be found by interpolation if needed. 6 There are mature algorithms for generating these meshes, and novel methods are still an active area of research—particularly regarding robustness in automated meshing [14, 15, 16, 17]. Ge- ometries derived from CAD models can be particularly problematic, as minor inconsistencies in the boundary representations lead to invalid tesselation of the volume. In addition to which polyhedra are used as elements, meshes are categorized by whether they are structured, and whether they are conforming. In conforming meshes there is no partial overlap between neighboring elements: there cannot be a node of one mid-edge or mid-face in the other. The cartoon mesh of figure 1.4 is thus nonconforming, as seen on the border of red and black re- gions. Since FEM’s derivation assumes a conforming mesh, most meshing algorithms are designed to never introduce nonconformities; satisfying this constraint will complicate mesh generation, as well as modification of existing meshes. Conversely, both structured and unstructured meshes see frequent use. Generally, an unstruc- tured mesh will be tetrahedral and a structured one hexahedral, though this is not definitionally true. Structured meshes specifically refer to those isomorphic to a regular grid, where the con- nection between elements is implied by their relative location within the computer’s memory. An unstructured mesh must explicitly track all connectivity between elements, nodes, and edges; it is thus more adaptable, but requires more computational resources to maintain and parse this topo- logical information. The mesh in figure 1.4 would seem to be unstructured, since the red elements do not lie on the same grid as the rest of the mesh; however, chapter 3 will demonstrate a represen- tation of such meshes that is neither truly structured or unstructured. 1.4.2 Conductivity discretization Once the mesh has been created, its topological nodes can be considered electrical nodes of a resistor network, with the discrete resistances calculated from the geometry of their containing elements. The admittance method is explicitly couched in these circuit-based terms[18]; though the finite element method is rarely presented as such, it too is mathematically equivalent to constructing a complex resistive network [19, 20, 21]. 7 1.4.3 Constraints & solving Any imposed voltages, or injected currents, compose the boundary conditions of the simulation. General FEM packages may perform fairly complex discretization of electrode surfaces, specif- ically constraining the resultant current density through each face; other approaches may simply treat the electrode as equipotential and merge all corresponding mesh nodes to a single electrical node. The boundary of the mesh is often set as ground to represent the field’s decay at infinity, and provide numerical stability if the simulation does not involve an explicit ground electrode. At this point, nodal analysis conveniently expresses the electrical network as a linear system of equations, which can be solved by typical linear algebra packages. With a voltage now calculated at each mesh node, potential in the interior of an element can be interpolated from its nodal voltages, for instance to determine the extracellular voltage at several neural compartments within a much larger element. 1.5 Potential improvements Presently, meshing and solving in the extracellular domain is far less optimized than the simulation of the intracellular domain. With neurons joined only at synapses, and compartments within a neuron linked only to the next and previous segment, it is relatively easy to divide an intracellular simulation among multiple threads, and even multiple machines [22, 23]. However, the same cannot be said for the extracellular space, where the meshing and solving processes are often slow, single-threaded procedures[24, 25]. As an approximation method, these mesh-based methods will generally rely on subsequent re- simulation using a finer mesh in order to verify that results are within a satisfactory margin of error. In a FEM approach, refining the mesh while maintaining conformity requires either that all ele- ments of the mesh are subdivided, or that significant modifications are made to transition smoothly from the newly refined elements to the much larger ones from the original mesh. While simpler to implement and safer to apply, global splitting multiplies the mesh’s element count and thus the size 8 of the system of equations. Though the conductance matrix is sparse, and thus practical algorithms for solving it are closer toO(N) than theO(N 3 ) of a general matrix inversion, this still results in a rapid expansion of the computational resources required as the mesh is refined. Further complicating matters, many neural simulations involve significant variation in their resulting field over the course of a simulation. Figure 1.5 illustrates a typical pattern of activity in a hippocampal slice following external stimulation, with drastically different distributions of spiking neurons at different timesteps. Creating a single mesh to adequately represent each step would require untenably high resolution, of which only a subset will be relevant at a given timestep. Figure 1.5: Simulated neural activity in a hippocampal slice following stimulation. Note the wide spatiotemporal variance over the course of the simulation. Reproduced from [26]. 9 Consequently, a means of quickly generating and modifying simulation meshes would help op- timize the use of computational resources in neural simulations. Early attempts will be discussed in the next chapter, focused on defining mesh manipulations in terms of equivalent circuits, an intuitive companion to the circuit theory employed in compartmental models. Though that spe- cific effort will prove untenable, it will set the stage for a more productive perspective on mesh generation that carries the remainder of the work. 10 Chapter 2 Equivalent-Circuit Transforms 2.1 Motivation Intuitively, the electrical properties of the extracellular space can be modeled as a 3-dimensional network of impedances. In practice, these will be simple resistors, since the reactive components of biological tissues are typically insignificant at the spatiotemporal scales relevant to neural activ- ity [27]. Later chapters will elaborate on the calculation of these discrete component values from the bulk material properties; for now, it is simply assumed that there is an existing resistor net- work, adequately representative of the real extracellular space, to be modified for computational efficiency. Of course, there must be a connection between the intracellular compartmental model and this extracellular network. The most obvious solution is to place a node in the 3-d network at the location of each compartment, electrically connected to the extracellular node of the compartmen- tal model. However, realistic models may contain millions of neurons, whose compartments are orders of magnitude smaller than the geometric features of the surrounding tissues [23]. Maintain- ing a mesh refined enough to capture these details, and large enough to encompass the simulation domain, introduces an extreme computational burden, necessitating compromises in the model’s detail. One potential approach is to adjust the mesh resolution per time step of the simulation, ac- cording to the present distribution of neural activity; quiescent segments’ nodes can be eliminated 11 to simplify the mesh, then re-introduced when required. To avoid introducing error (or excessive calculation that outweighs the time saved), the operation should be local; that is, the solution in the remainder of the mesh should be independent of the topology changes in the area under refinement, since the bulk conductivity is physically independent of activity. This chapter seeks to derive this relationship between coarsened and refined meshes in terms of formal circuit theory, a conceptual framework already familiar to most engineers. 2.2 Method 2.2.1 Setup For simplicity, we will begin with a homogenous, 2-dimensional square region as the element to be modified, rather than proceeding to the more applicable case of a 3-dimensional object. By symmetry, a square element has the same equivalent resistance Z edge as measured between any two adjacent vertices, and the same Z diag between nodes on the diagonal. Figure 2.1 indi- cates these quantities on a standard electrical schematic. In the circuit theory sense, these vertices represent the terminals of the resistor network; if the network is altered by adding, removing, or changing the value of the individual resistances, but the net resistance between terminals remains the same, the two variants of the network are said to be electrically equivalent (or are equivalent circuits). Figure 2.2 illustrates one potential subdivision of a square element, whoseZ diag ought to remain equal. 2.2.2 Solver From elementary circuit theory, the elementary equivalent-circuit operations R ′ = ∑ R i (series) R ′ =( ∑ R − 1 i ) − 1 (parallel) 12 Ω Ω Figure 2.1: Characteristic impedances of a homogenous 2d square Ω Ω Figure 2.2: A region before and after refinement. Resistor values post-refinement should be chosen such that net impedance stays the same. R ij =R i0 R j0 N ∑ k 1 R k (Y-∆) alter the resistor network but maintain its equivalent resistances. For networks whose topology is a planar graph (i.e., can be drawn without overlapping edges), it is guaranteed that some combi- nation of these transforms will reduce the network to a single equivalent resistance between any two nodes [28]. While possible to carry out these transformations by hand, the more compli- cated 3-dimensional networks will prove tedious and error-prone. Thus, an automated approach is preferable. A custom netlist solver is implemented in MATLAB with the Symbolic Math Toolbox; all novel code is available at https://github.com/cbcGirard/SymbolicNetSolver. The circuit in question is represented in typical netlist format, with every electrical junction given a numeric 13 tag, and each resistor specified by its value and the tags of its terminals. After specifying the pair of terminals under test, the solver attempts to eliminate the other nodes of the circuit by repeated applications of the series/parallel and Y-∆ transforms, as outlined in 1. The solver will then return a closed-form, symbolic expression for the net resistance between the terminals in terms of the original network’s resistors. Algorithm 1 Equivalent resistance solver ⊕ :append ⊖ :delete A\B :subsetof AnotinB procedure REDUCENETWORK(R list ,nodes,terminals) while LENGTH(nodes)> 2 do fornode∈{nodes\terminals} do R ′ ← r inR list |node∈r.nodes if LENGTH(R ′ ) == 3 then R list ⊕ YDELTA(R ′ ) ▷ Y-∆ R list ⊖ R ′ nodes⊖ node else if LENGTH(R ′ ) == 2 then R list ⊕ ΣR ′ ▷ series R list ⊖ R ′ nodes⊖ node else fornode ′ ∈{R ′ .nodes\node} do R ′′ ← r inR ′ if{r.nodes\node}==node ′ ▷ parallel R R list ⊕ 1/(Σ(1/R ′′ )) R list ⊖ R ′ end for end if end for end while end procedure 14 2.3 Results 2.3.1 2d square: electrical equivalence is impossible The values of resistance for pre- and post-split networks are illustrated in figure 2.2. The two characteristic impedances of the refined network are determined by the automated solver to be Z diag,1 =k E R 0 2k C R 0 +k E R 0 k C R 0 +k E R 0 (2.1) Z edge,1 = k E R 0 2 3k C R 0 + 2k E R 0 k C R 0 +k E R 0 (2.2) while their pre-refinement equivalents, Z diag,0 = R 0 and Z edge,0 = 3 4 R 0 respectively, can be de- termined by inspection. Setting the pre- and post-refinement values equal yields the system of polynomial equations k 2 E − k E + 2k E k C − k C = 0 (2.3) 4k 2 E − 3k E + 6k E k C − 3k C = 0 (2.4) which are only simultaneously valid for the trivial solutions (k E = 0.5,k C =∞) or (k E = k C = 0). Figure 2.3 plots the equations for equal Z diag and Z edge independently, illustrating that the two cannot be simultaneously satisfied except at the the trivial solution (k E = 1 2 ,k C =∞), which is degenerate to the original, un-refined network. The other solution is, of course, physically meaningless. Thus, there is no set of resistor values that will make the two networks electrically equivalent. 15 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 k E -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 log 10 k C edge diag Figure 2.3: k C as a function of k E , as would maintain equal impedance across the edge (blue) or diagonal (red) 2.3.2 2d square: error minimization Though true equivalence is impossible, a refinement strategy that minimizes error may instead be suitable. Denoting the unrefined network by subscript 0, we formally define the looser metric as the sum of errors Error= ∑ (Z 0 − Z 1 ) 2 (2.5) =(Z diag,0 − Z diag,1 ) 2 +(Z edge,0 − Z edge,1 ) 2 (2.6) Inserting the values from 2.1 and 2.2 expands to the rather unwieldy 100k 2 C k 2 E − 100k 2 C k E + 25k 2 C + 112k C k 3 E − 156k C k 2 E + 50k C k E + 32k 4 E − 56k 3 E + 25k 2 E 16(k C +k E ) 2 (2.7) This SSE is plotted below (wide view for values 1e-3 to 1e3 at left, and detailed view of the minimum at right). The error is relatively insensitive to k C , with 0.5<k E < 1 giving the lowest error. 16 (a) Wide view of log 10 Error vsk E andk E (b) Detail around minimum of error Figure 2.4: Refinement error in split mesh by resistor values Alternatively, expressing the SSE as 2 56k 2 E − 78k E + 25 k E k C + 25(1− 2k E ) 2 k 2 C + 32k 2 E − 56k E + 25 k 2 E 16(k C +k E ) 2 suggests potential minima, since the polynomial terms (56k 2 E − 78k E + 25) (1− 2k E ) (32k 2 E − 56k E + 25) have the respective roots k E =( 1 2 , 25 28 ) k E =( 1 2 ) k E =( 7± i 8 ) Ignoring the complex roots, candidate minima for error occur at SSE(k E = 0.5)= 5 64 k C + 1 2 2 SSE(k E = 25 28 )= 3025k C 2 196 + 15625 38416 16 k C + 25 28 2 17 As seen in figure 2.5, the former equation suggests error approaches zero as k C →∞; this is degenerate to the original mesh. The latter (withk E = 25 28 ≈ 0.893) has a fairly steady error around 3% fork C < 0.1 (though its minimum occurs atk C = 25/847≈ .0295). 2.3.3 Nodal interpolation without splitting Since the previously discussed splitting method is incapable of producing an electrically equivalent result, an alternative circuit topology is proposed. Illustrated in figure 2.6, an interior node can be introduced without splitting any of the extant edges, but rather by splitting the square into triangular subsets. Analysis of the resulting networks proceeds as in the previous section. This time, the systemZ diag,1 − Z diag,0 = 0,Z edge,1 − Z edge,0 = 0 expands to 2k C k E − 2k C − k E = 0 (2.8) 8k C k E (3k C +k E )− 18k C k E − 24k C − 3k E = 0 (2.9) with trivial roots(k C = 0,k E = 0) and(k E = 1,k C =∞). As in the previous analysis of the refined mesh, relative error in impedance can be quantified by the sum of squared errors SSE =    1600(k E − 1) 2 k 4 C + 32k E 28k 2 E − 103k E + 75 k 3 C +4k 2 E 32k 2 E − 268k E + 325 k 2 C − 4k 3 E (28k E − 75)k C + 25k 4 E    16(2k C +k E ) 2 (4k C +k E ) 2 (2.10) From the SSE, the polynomial terms (k E − 1) (28k 2 E − 103k E + 75) (32k 2 E − 268k E + 325) (28k E − 75) 18 10 -4 10 -2 10 0 10 2 10 4 k C 10 -10 10 -8 10 -6 10 -4 10 -2 10 0 SSE k E =0.5 k E =25/28 (a) Error at potential 10 -4 10 -3 10 -2 10 -1 k C 0.03 0.032 0.034 0.036 0.038 0.04 SSE (b) Detail near minimum Figure 2.5: Refinement error at candidate minima of split mesh 19 Ω Ω Figure 2.6: Schematic of diagonal refinement with roots k E = 1 k E =(1, 75 28 ) k E = 67± √ 1889 16 k E = 75 28 suggest possible minima at SSE(k E = 1)= 356k C 2 + 188k C + 25 16 8k C 2 + 6k C + 1 2 SSE(k E = 75/28)= 220900k C 4 49 − 11250000k C 2 2401 + 791015625 614656 16 8k C 2 + 225k C 14 + 5625 784 2 Again, error only approaches zero in the trivial case(k E = 1,k C =∞) degenerate to the original mesh. Otherwise the smallest error (≈ 0.155%) occurs for (k E =(67− √ 1889)/16≈ 1.471,k C ≈ 1.458). 2.4 Discussion and alternative formulations While intended to be a first step in deriving a theoretical basis for dynamically adding nodes to a mesh, the present study instead suggests that altering the topology within a resistor network 20 (a) log 10 error vsk C andk E (b) Detail of error near potential minima 10 -4 10 -2 10 0 10 2 10 4 k C 10 -10 10 -8 10 -6 10 -4 10 -2 10 0 10 2 SSE k E =1 k E =75/28 k E =(67+1889 0.5 )/16 k E =(67-1889 0.5 )/16 (c) Error as function ofk C at candidate minima. Figure 2.7: Refinement error in crossed mesh by resistor values 21 will alter its equivalent impedances, regardless of how the individual resistances are changed to compensate. In the simplest cases, it is possible to derive the conditions of equivalent impedance across individual pairs of nodes, but each impedance presents mutually incompatible constraints. There are technically optimal solutions, in the sense of minimizing the sum of errors in these impedances, but the minima are shallow and seemingly divorced from physical significance. 3-dimensional cases will pose additional challenges, as their non-planar structure may prevent even the derivation of a single equivalent impedance by the same method employed here; though valid star-mesh transforms may be possible in these cases, they are not guaranteed to exist as in the 2-d cases investigated here [28, 29, 30]. Further complications such as anisotropy, or other phenomena expected in a practical application, are unlikely to help in arriving at a solution. Though it would be of theoretical interest to confirm that equivalence is not possible in those cases, such a negative result would be of little practical use. Instead, this line of inquiry is suspended in favor of a single approximation rule, applied uni- formly regardless of whether a cuboid represents an ”original” element of the mesh or the result of later subdivision. Three candidate formulations are illustrated in figure 2.8 and described in the following sections; their relative strengths and weaknesses will be empirically determined in the next chapter. 22 Standard Admittance Graph-Dual Admittance Trilinear Finite-Element Figure 2.8: Circuits resulting from different methods of discretizing conductivity in a cuboid region 23 2.4.1 Admittance method The most intuitive the admittance method introduced in [18]. In short, the orthogonal components of current flow are taken to be uniform between opposing faces, passing through a conductance of σ A l . This net conductance of the cuboid is divided evenly among its four parallel edges, equivalent to the resistor network shown in figure 2.9. By redistributing conductance to the edges in this man- ner, the connectivity of the electrical network is thus given by the already-established connectivity of the mesh. Because anisotropic conductivity takes the form σ =       σ x 0 0 0 σ y 0 0 0 σ z       (2.11) handling anisotropy within the admittance element simply requires using the proper component of the tensor for each edge orientation. Figure 2.9: Schematic of admittance method voxel In contrast to the more common finite element method, this approach does not intrinsically require that the mesh remain conforming. Allowing nonconforming elements allows splitting only a subset of elements without needing to ”repair” the resulting irregularities, greatly simplifying the 24 refinement process compared to traditional algorithms. Unfortunately, the literature appears silent on how these irregularities may affect the accuracy and numerical stability of the subsequent simu- lation. Since the present work seeks a fast remeshing process, it will employ this simplification of the nonconforming mesh; as a corollary, the behavior at nonconforming nodes will be of particular interest, in case they exhibit particular instability. 2.4.2 Mesh dual Alternatively, for the same net conductance of σ YZ X between the faces normal to ˆ x, instead of decomposing it into σ YZ 4X along each edge of the element, the faces can remain as the electrical nodes of the system; splitting to add a node at the center of the element yields two conductances of 2σ YZ X , and likewise for the other two axes. This configuration is shown in figure 2.10a, and the topology of the resulting network across elements is effectively the graph dual of the original mesh. Choosing this representation instead of the classic admittance method allows a natural solution to mesh nonconformities: either the smaller elements’ nodes are merged to the larger element’s, as shown in figure 2.10b, or else the larger element’s conductance is split among the smaller ones, as shown in figure 2.10c. Either one enforces continuity of current between adjacent elements as a consequence of the electrical topology; in contrast, the edge-only conductance of the admit- tance method must shunt current laterally when it is incident on a face. Though the merging is algorithmically easier (simply replacing nodal tags), it results in a loss of spatial resolution. Thus, the latter formulation of splitting will be compared as the nearest equivalent to the traditional ad- mittance method. This approach does, however, require an extra step of determining adjacency between elements of different sizes, since their facial nodes are not geometrically coincident. A similar topology was shown to form a lower bound on conductivity compared to the upper bound of the trilinear finite element, at least in the scenario of a uniform grid [19]. This suggests the possibility of obtaining a reliable confidence interval on the true solution from a single mesh, 25 (a) Equivalent circuit of an element’s dual (b) Merging nodes at nonconforming junction (c) Splitting conductance at nonconforming junction Figure 2.10: Equivalent circuits of the mesh dual 26 by sequentially employing the two discretizations. If so, it could be used to eliminate a additional remeshing step to verify convergence, avoiding a much costlier operation. 2.4.3 Trilinear finite element Typical presentations of the finite element method are couched in terms of nonspecific differential equations, and would seem entirely foreign to the resistor network concept underlying the admit- tance method. Though this flexibility can be beneficial, it also imposes significant complexity in implementation. However, reinterpreting its results specifically in terms of circuit theory suggests a relatively small computational overhead, while its theoretical basis implies a greater ability for the element to capture the behavior of the potential field. For electrical conductance, the governing differential equation is Laplace’s equation 1 ∆φ = 0 (2.12) which the finite element method seeks to approximate on the discrete elements. Full derivation of the general case is outlined in appendix B, but the net result is an expression for conductance between nodesi and j g i,j = Z Ω σ∇h i · ∇h j dA (2.13) where h n functions are basis functions on the vertices, whose sum provides interpolation of any point within the element in terms of potential at the vertices. In the case of a hexahedral element, trilinear interpolation can be considered the sum of basis functions and nodal voltagesv i u(x,y,z)= 7 ∑ i=0 v i h i (x,y,z) (2.14) where bothu andh i are polynomials of the form a+bx+cy+dz+exy+ fxz+gyz+hxyz (2.15) 1 Essentially the continuous case to the discrete Kirchoff Current Law, since∆v=∇· σ∇v= divi 27 The coefficients of the h i terms can be determined by solving the system          1 x 0 y 0 z 0 x 0 y 0 x 0 z 0 y 0 z 0 x 0 y 0 z 0 1 x 1 .............................. ..................................... 1 .......................... x 7 y 7 z 7                a i . . . h i       = b (2.16) where b n =        1 i= j 0 i̸= j (2.17) By symmetry, conductances within the element fall into one of three categories: along an edge (one differing coordinate between the nodes) across a face diagonal (two differing coordinates) across the element (three differing coordinates) If the element’s sides are orthogonal, the integral in 2.13 can then be evaluated analytically, which yields: g x = YZ 9X σ x − XZ 18Y σ y − XY 18Z σ z (2.18) g xy = YZ 18X σ x + XZ 18Y σ y − XY 36Z σ z (2.19) g xyz = YZ 36X σ x + XZ 36Y σ y + XY 36Z σ z (2.20) For an isotropic cube, the latter two reduce to ℓ 0 12 and the first to zero. For a mesh of uniform cubes, these elements theoretically provide a closer approximation of the underlying conductivity than the equivalent admittance elements[19]. Moreover, the trilinear 28 Figure 2.11: Equivalent circuit of a trilinear finite element interpolation underlying the finite element effectively models variation in current density within the element, since its gradient is itself a quadratic function. In contrast, the admittance element assumes a constant current vector over its volume in order to separate its spatial components; in this regard, it is similar to the typical tetrahedral elements known to underperform their hexahedral counterparts [17]. 2.5 Conclusion Though conceptually, subdividing space to attain a finer resolution leaves its electrical properties unaffected, in practice this will intrinsically change its discrete representation. No set of resistance values in the modified region can achieve electrical equivalence with the original, and so any means of embedding a new node will alter the effective circuit. Though it is thus impossible to implement an equivalent-circuit transform as a distinct oper- ation, the adaptation process can instead be understood as a re-application of the same logic that created the original circuit. The more complicated finite-element formulation is theoretically more accurate per-element than the simpler admittance method and its graph-dual formulation, but it is unclear whether this advantage extends to nonconforming meshes that violate its underlying as- sumption of continuity. Similarly, theory suggests that graph-dual formulations provide a lower 29 bound on conductance compared to those based on the native mesh topology, but only under the assumption of continuity. Consequently, the following section will take a different tack: rather than seeking a distinct operation to quickly insert detail in an existing mesh, it will present an organizational scheme for rapidly generating and adjusting the mesh itself. 30 Chapter 3 Adaptive Octree Meshing 3.1 Motivations Further investigation in the literature suggests it is unusual to consider adaptation as a different set of mathematical operations from those employed in generating a simulation mesh de novo. Rather, the geometric element is taken as the basic unit, and its effective resistances are determined solely by its geometry. In this framing, adaptation seeks to adjust the mesh geometry so that its discretization rule better approximates the true solution, generally by shrinking the elements’ size where the solution changes most rapidly. Efficient mesh adaptation, then, is intertwined with efficient mesh generation. The actual data structures that encode the mesh geometry must allow easy querying of specific regions; insertion, deletion, and modification of the elements within the region; and yet maintain coherence with the rest of the mesh. In modern computing, parallelism is key to increasing performance; even low-end processors are generally multicore, and at the high end, it is far easier to add nodes to a cluster than to boost raw speed by a similar percentage. However, few computations are trivially parallelized, including those involved in meshing: local changes in topology such as splitting an element affect neighbor- ing elements, a minor concern if said neighbors are easily identified, and potentially unworkable if the entire mesh must be searched element-by-element. 31 Using subdivision of a cuboid as a fundamental operation of both refinement and adaptation, we can arrive at a mesh structure that satisfies the needs outlined above. This octree structure enables modification of the mesh by independent threads, spreading the computational burden, while these local operations are inherently coherent with the global mesh because it mathematically encodes their location. 3.2 Methods 3.2.1 Overview Figure 3.1 gives a conceptual overview of the simulation process. Starting with the entire simula- tion domain as an element, elements are recursively split if they are larger than a target size that is a function of their location. Once all elements satisfy this size criterion (or after a maximum number of iterations is reached), their local conductances are collected into a global netlist. If multiple nodes of the mesh lie within an equipotential region, such as an electrode, they are merged into a single point. Finally, this effective netlist is converted to a linear system of equations and solved for the potential at each node, based on any currents injected into the mesh and applied voltage constraints. 3.2.2 Octree-based meshing This top-down, recursive subdivision results in a structure known as an octree. This differs from the unstructured meshes typical in finite-element applications, where the order of nodes and elements is essentially arbitrary and connectivity must be explicitly recorded. However, it is also distinct from so-called structured meshes, where the mesh is essentially a regular grid mapped to the sim- ulation geometry, and thus arrays implicitly reflect the connectivity. Similar mesh arrangements have been previously used in neural simulations, though generated instead through a bottom-up process of merging voxels [18, 31]. In those studies, it appears that the mesh structure is not in- voked after construction, instead converting each resulting hex element to an equivalent resistor 32 ... Figure 3.1: Constructing a simple mesh. 33 network; thus, it would be considered a type of unstructured mesh under the formal definition, even if the resulting arrangement is similar to that proposed here. The octree structure is particularly amenable to binary representation. Each octant can be mapped to an integer between 0 and 7 according to its position in Cartesian coordinates, with each axis represented by one bit as shown in figure 3.2a. A list of such integers will thus uniquely define any cubic sub-region as illustrated in figure 3.2b; when an element is split, each child’s global identity is the integer tag of its octant appended to its parent’s list of integers. 6b110 4b100 0b000 1b001 3b011 7b111 5b101 2b010 (a) Octant numbering by Cartesian position 1 1,4 1,4,7 (b) Element identification by list of tags Figure 3.2: Diagram of octree structure 34 Similar logic applies to identifying elements’ vertices, except after n rounds of subdivision there will be 2 n elements along an axis but 2 n + 1 nodes. Consequently, a numerical tag can be derived for each vertex from its position within the element, and the identity of the element itself. Again, this number will be unique to that point in the simulation domain, so that adjacent elements will return the same integer for a shared node. Since this global identifier is in terms of local quantities, element-specific calculations (such as discretizing conductivity) can happen in parallel, without the need to later reconcile individual node numberings; moreover, because these operations are all integer-based, this paradigm avoids floating-point errors that could mis-tag a node and break mesh connectivity. 3.2.3 Splitting metric As a rule of thumb, accuracy is most sensitive to mesh resolution in regions where the gradient of the solution is highest. In practice, this relationship is not rigorously quantified; each simulation is a unique combination of geometry and boundary conditions, whose solution is by definition not knownapriori. Common methods of inferring a simulation’s reliability employ re-simulation using successively finer meshes to empirically estimate the solution’s convergence. Reliability of these approaches generally come at the cost of computational inefficiencies. Es- pecially in traditional meshing pipelines, where mesh conformity is required, it is often simpler to subdivide every element in the mesh for the next approximation, increasing the required com- putational resources exponentially with respect to the number of iterations. Otherwise, significant effort is required to smoothly transition the size of elements between the refined region and the rest of the mesh, at which point an entirely new mesh may be more practical. Thus, having a heuristic for element size is still useful at the point of initial mesh generation. Generally, there will be some function of location within the mesh f(x,y,z) that returns the char- acteristic lengthℓ 0 , a measure of the desired element size at that location 1 . This size field can be 1 Exactly what geometric quantityℓ 0 refers to seems to vary by author, aside from the fact that it is a 1-dimensional length. In this work, it is defined as the geometric mean of the cuboid’s edge lengths; regardless, its utility is as a relative measure of size, not necessarily an absolute one 35 an arbitrary function or set of functions, where the smallestℓ 0 is used to define the element at that point [32]. For electrophysiology studies, simply using a scalar multiple of distance from a source is a practical choice, due to the 1/r dependence of the analytical solution for a point or spherical source [18]. In this work, the size field takes the form ℓ 0 =α∥r∥=∥r∥2 − kN (3.1) where∥r∥ is the Euclidean distance from the element’s center to the source’s center,N is the maxi- mum number of subdivisions permitted, andk controls mesh density, ranging from fully subdivided everywhere (k= 1) to as coarse as possible while attaining a full split at the source (k= 0). 2 Unless otherwise noted, a default value ofk= 0.2 is used throughout, as it was empirically determined that lower values degrade accuracy, while higher values increase computational burden with relatively small increases in accuracy. 3.2.4 Netlist creation, boundary conditions, and solving The present work essentially follows prior literature on the admittance method [18]. Some of the stepwise logic presented in the literature is instead implemented here as sparse matrix operations, taking advantage of extant libraries optimized for parallel execution. In short, each leaf 3 contributes the 12 discrete conductances between its nodes to the global network, resulting in a pair of arrays: [g] conductances between the pairs of nodes [i, j]. This is equivalent to a sparseN 0 byN 0 sparse matrixG 0 in coordinate list format, since the net conductance between nodesi and j is simply the sum ofevery discrete conductance between those nodes. 2 A full attempt to derive the ”true” function that maps(0,1) from minimally to maximally refined meshes can be found in A. The mathematical expressions are too cumbersome to derive a closed-form solution, but it is at least shown that equation 3.1 at least produces the desired meshes at the extreme values of k, even if more moderate values result in the same mesh. 3 I.e., terminal octant that has not been split 36 At this point, equipotential nodes can be merged to the same node tag, reducing the size of the matrix to N 1 by N 1 . Thus, for an electrode injecting current into the extracellular space, the conductive network formed from the mesh geometry encodes how the current distributes through the mesh, but the solving process only needs the total current through the electrode. While these node indices are being rearranged, it is convenient to group those with floating potential, with injected current, and fixed voltage in that order, to simplify later operations. Rearranging Kirchoff’s Current Law at node n, with any current injected into the node from outside the mesh asi ′ n , results in ∑ i n = 0 (3.2) N 0 ∑ j=0,̸=n g nj v j − v n +i ′ n = (3.3) v n N 0 ∑ j=0 g nj − N 0 ∑ j=0,̸=n g nj v j =− i ′ n (3.4) G 1 v= i (3.5) where G ′ 0 = G 0 + G ⊤ 0 (3.6) G 1 = diag G ′ 0 [1,1,...,1] ⊤ − G ′ 0 (3.7) Nodes with known voltages (such as grounded periphery or voltage-controlled electrodes) can be eliminated from the system by defining v 0 =[0,0,...,v N 1 − 1 ,v N 1 ] and performing the substitution i ′ = i− G 1 v 0 (3.8) 37 truncated to the first n=N 1 − n v entries (the floating and/or current-injected nodes with unknown potentials). Similarly trimming G= G 1 (0 :n,0 :n) yields the final linear system of equations Gv= i ′ (3.9) Due to the extreme sparsity of the conductance matrix (generally 7 nonzero entries per row), iterative methods are significantly faster than explicit inversion of the matrix in all but the coarsest meshes. In this study, the conjugate-gradient method is used, as described in prior literature on the admittance method [33, 18]. Other Krylov subspace methods have been reported in similar applications [34]. 3.2.5 Validation tests 3.2.5.1 Common features All simulations are executed on a consumer desktop machine (Ryzen 9 3900x, 32 GB RAM) running Kubuntu 20.04. As a balance between performance, ease of development, and flexi- ble reuse, the code is developed as a Python package incorporating several open-source libraries [35, 36, 37, 38, 39, 40, 41, 33]. Future release is planned at https://github.com/cbcGirard/ xcell. In each simulation, the source is a 1µm sphere at which sufficient current is injected to theo- retically produce 1.0V potential at its surface; in other words, I src = 1.0V∗ (4πσr src ) (3.10) The simulation domain is a cube 200µm per side, with homogenous isotropic conductivity of 1.0S/m 4 . 4 Though chosen for numerical convenience, this is roughly the correct magnitude for physiologically relevant tissues [42] 38 At this size of simulation domain, the theoretical potential at the boundary is r elec /r node ≤ 1%. These theoretical values can be easily set as the boundary conditions to mitigate error due to boundary effects; however, in practice these effects are so minor that simply grounding the nodes suffices. To benchmark the performance of the octree mesh, at each stage of refinement the simulation is repeated using the next-largest (in terms of element count) uniform grid. 3.2.5.2 Conductance discretizations The admittance, FEM, and mesh-dual approaches described in 2.4 can be easily substituted for each other in the course of a simulation. Thus, after constructing a set of octree meshes at various resolutions, the relative performance of each discretization can be directly compared. Because the mesh dual’s electrical nodes no longer correspond to mesh vertices, there is no explicit expression for potential at the center of the domain. To provide a closer equivalent to the other meshes, the potential at the vertices of each element is additionally interpolated from the values at the dual nodes according to u(x,y,z)=ax+by+cz+dx 2 +ey 2 + fz 2 +g (3.11) with coefficients derived from                 − 1 0 0 1 0 0 1 1 0 0 1 0 0 1 0 − 1 0 0 1 0 1 ....................... 0 0 1 0 0 1 1 0 0 0 0 0 0 1                                 a b c ... f g                 =                 v − x v x v − y ... v z v 0                 (3.12) 39 Since coarser meshes generate no dual nodes that lie within the source, unlike the center vertex present in all meshes, the site of current injection is instead moved to the closest dual node. 3.2.5.3 Error metrics Defining a metric for error equally applicable across domain sizes, mesh types, and element reso- lutions is surprisingly non-trivial. Error at a node is easily found by the difference between theo- retical and simulated potentials at that point; however, simple methods of aggregating these errors (such as root-mean-square) assume nodes are all equally valid. Though somewhat reasonable for a uniform grid, where it can be considered a uniform sampling of the actual field, it will not suffice for the varying resolution of the octree mesh. Size of the simulation domain also complicates matters. The total number of nodes on the periphery will tend to increase with distance, and these nodes will exhibit the smallest potentials (and thus errors). Attempting to normalize an error metric by the number of nodes (as in RMS) will thus artificially deflate the metric by overrepresentation of these nodes; on the other hand, using relative error is likely to inflate it due to the small divisor of the theoretical values. As a compromise between these concerns, the following is proposed as a meaningful quan- tification of a simulation’s error. The analytical solution at distance r from the source’s center is v analytic (r)=        1/r r>r src 1 r≤ r src (3.13) Error is calculated at each node, and used to construct a piecewise linear function Error(r). For a simulation domain whose furthest point isx max from the origin, net error is given by E sim = Z x max 0 |Error(r)|dr (3.14) 40 conveniently calculated by trapezoidal integration. This may be normalized by V analytical = Z x max 0 v analytical (r) dr (3.15) =v src r src + ln x max r src (3.16) to compare simulations with different source amplitudes, and to give a more intuitive sense of error as a fraction of the true solution. 3.3 Results 3.3.1 Discretization formula Results from a representative sample of meshes are shown in figure 3.3; log-log axes are used to better visualize how each discretization’s potential field varies around the true analytical solution. The shift in effective source location for duals of coarse meshes result in wide variance among points of the samer, depending on which direction the actual source moved; the FEM simulations show variance concentrated near nonconforming nodes, as also seen in the admittance method. Figure 3.4 shows the net accuracy of the three discretization methods for the same mesh. The computational requirements of each are shown in figure 3.5. The finite element approach requires moderately more time in discretization and solving than the admittance method does; calculation of the mesh dual is a non-trivial contributor to simulation time, but the increase in solver time has a still greater impact. As the fastest and best-performing formulation, the traditional admittance method will be used for subsequent investigations. 3.3.2 Octree and uniform meshes Figure 3.6 shows a planar slice of a simulation using a moderately-refined octree mesh, illustrat- ing a typical structure. Note that at this stage, multiple mesh nodes lie inside the boundary of the source; these are constrained to be equipotential, and thus behave as a single electrical node. 41 However, this degree of refinement does result in explicit representation of the source boundary in the mesh, unlike coarser cases where a single node implicitly represents the entire source. Quantitatively comparing such images is difficult, so figure 3.7 presents simulated solutions and error as a function of distance from the origin (and thus, source). In each row, the uniform mesh has approximately the same number of elements as the octree mesh, whose resolution is governed by the maximum recursion depth allowed. Note that, as shown in the bottom plots of each row, all uniform meshes are so coarse that the source is represented by a single node, but the latter two octree meshes have achieved explicit representation of the surface. Additional detail of a typical error distribution can be seen in figure 3.8. 42 Figure 3.3: Solutions from admittance, FEM, and mesh dual methods. The mesh dual is qualita- tively ”smoother” than the other methods, but does not consistently form a lower bound on the true solution as anticipated. 43 Figure 3.4: Accuracy of admittance, trilinear FEM, and dual methods applied to the same octree mesh 44 Figure 3.5: Total computation time of simulation as a function of mesh size. Though achieving similar or better accuracy, the standard admittance method is the simplest and thus fastest to cal- culate. 45 Figure 3.6: Solution and error in XY plane (Z=0) for octree depth 10 46 Figure 3.7: Error distribution of uniform (left) and octree (right) meshes at maximum depths of 5, 9, and 17. Top traces of each set show simulated and analytical voltages, with the absolute error between them in the center trace. The bottom trace shows the corresponding distribution of element sizes. 47 Figure 3.8: Detail of error at non-source nodes in an octree mesh, maximum depth of 14 48 Finally, the net error of the uniform and octree meshes are shown in figure 3.9, where the octree mesh consistently outperforms the similarly-sized uniform mesh. Figure 3.10 compares the computational requirements of each step in the simulation process. Wall times appear linear with respect to element counts, though CPU times increase more rapidly as parallelized routines occupy a growing portion of the simulation. The actual generation of the mesh geometry is the slowest step for both mesh types, and takes several times longer in the octree mesh than the uniform. Figure 3.9: Net error of uniform and octree meshes by element count (top) and smallest mesh element (bottom). Note that net accuracy is better predicted by resolution near the source, though the octree mesh will use far fewer elements to achieve the same central resolution. 49 (a) Contribution of each simulation step to total execution time. All exhibit O(n) complexity, so linear extrapolation to larger simulations is reasonable. (b) Total simulation time by peak mesh resolution. Though uniform meshes are faster per-element to con- struct as seen above, octree meshes are orders of magnitude faster in achieving a given spatial resolution. Figure 3.10: Time requirements of simulations 50 A separate set of simulations, varying in source radius, initial domain size, and mesh density, are aggregated in figure 3.11. By perturbing size of the domain and the meshing parameters, each tracks the changing behavior as mesh resolution approaches source radius, after which the source is represented by multiple nodes. All follow the same basic trajectory: error decreases monotonically untilℓ 0 /r src =π. At this point, the simulated source voltage (which has steadily increased) agrees with the theoretical value; it then overshoots dramatically until re-converging as the mesh better approximates the surface of the source. Figure 3.11: Overshoot as localℓ 0 approaches source radius. Trace colors represent combinations of domain sizes (50 and 500µm), source radii (1 and 5µm), and mesh density factors (0.2 and 0.4), which all follow the same trajectory.ℓ 0 /r src =π is indicated by the dashed line. 51 3.4 Discussion In short, the octree mesh outperforms the uniform grid in all tested cases, as seen in figure 3.9; the largest meshes approach the upper limit of the test machine’s resources, at which point the uniform mesh still represents the source as a solitary node (cf. bottom row of figure 3.7). 3.4.1 Singularity effects The fact that error in the octree mesh is minimized atℓ 0 /r src =π (figure 3.11) is of limited practical use as a refinement target, due to the low likelihood of a realistic simulation achieving those specific dimensions. It does, however, suggest that error, at least as quantified here, is dominated by regions closest to the source, and that care should be taken when elements’ℓ 0 are slightly larger than the source’s radius; further subdivision to introduce multiple nodes within the source boundaries is likely warranted in this case. 3.4.2 Nonconforming elements In the octree meshes, error appears most prominently at the boundaries between elements of dif- ferent size, as illustrated in figure 3.8. These points correspond to non-conforming portions of the mesh, where the regularity of connections between the nodes is disrupted. Little has been pub- lished about the consequences of mesh nonconformance per se, aside from empirical comparison with a conforming mesh[43, 44, 45]. Here, it appears that the nonconformity indeed introduces error into the solution, though offset by the reductions won by the increased resolution within. Thus, these nonconforming points suggest a target for improving the method. In prior studies using the admittance method, element sizes within a mesh tend to span only a few subdivisions [46, 47, 48, 31]. Use of a top-down meshing methodology instead could result in far deeper levels of splitting, and thus stronger motivation to smooth the effects of these discontinuities. This possibility is explored further in the next chapter. 52 3.4.3 Simplifications and shortcomings As presented, there is no mechanism to prevent splitting of elements already fully contained by the source, whose child vertices will ultimately be merged into the single source node anyway. This is reflected in the relative plateauing of error in figure 3.9 at higher element counts, as the additional elements do not contribute to increased accuracy. Though this over-refinement does not affect the solver time, it does overstate the computational burden of the other steps. Such deep levels of refinement (with thousands of elements inside the source) are not likely to arise in practical simulations, so such a modification would be more a safeguard against unreasonable parameters than a genuine optimization. Similarly, each mesh is generated de novo in order to fairly represent the relative performance between the uniform and octree structures per se. However, the motivation of the octree structure is its ease of adaptability; in practical use, the computational cost of adding resolution to an extant mesh will be a fraction of the cost to generate the same mesh anew, but changing the uniform resolution requires a new mesh. Moreover, the simulations as implemented favor the uniform mesh; its process relies on optimized and compiled library code, while the novel octree code has additional overhead due to the Python interpreter 5 . As a future goal, re-implementation of the meshing logic in a compiled, optimized program is expected to greatly accelerate the octree method. The method also makes several simplifying assumptions about electrode behavior. In reality, actual electrode-electrolyte interfaces involve significant capacitive coupling by design, avoiding charge transfer between electrode and solution that could lead to toxic or degradative chemical reactions[49, 50]. Combined with the technically nonzero resistance of the electrode material, treating the electrode surface as a fixed potential is necessarily an approximation[11]. However, it is expected that few scenarios warrant the added complexity of explicitly mod- elling these characteristics in simulating the extracellular space itself. Practical studies suggest 5 As implemented, the performance gap is lessened somewhat by the Numba library’s just-in-time compilation of some functions[40]. Unfortunately, recursive functions and self-referential objects are not supported at this time, but are central features of the octree generation. 53 better numerical stability when the range of finite conductivities remains within a few orders of magnitude, a boundary stretched by explicitly including both electrode and physiological values [51]. If desired, the capacitive effects can be captured by a lumped model such as the Randles circuit or variants thereof, treating capacitance separately from the purely-resistive extracellular space[52, 50]. Finally, the same meshing algorithm can be used with complex admittances in- stead of pure conductances if required, at the cost of a more complicated discretization and solving process. Finally, the scope of this study is limited to meshes of hexahedral elements, though tetrahedral or mixed hex-tet approaches remain more common. It is already well-established that hexahedral meshes are more accurate on a per-element basis[26, 17]. Tetrahedral approaches are still used in part because of their maturity; robust tetrahedral algorithms exist for nearly any geometry, while automated generation of conforming hexahedra is still an active research topic[14, 53, 15, 16]. If the requirement for mesh conformity can be relaxed, the octree approach thus circumvents the issue of conformity to deliver a usable hexahedral mesh. 3.5 Conclusion Constructing a simulation mesh by recursively splitting the simulation domain results in a structure well-suited to later manipulation. Though somewhat slower to construct than a uniform grid, it provides greater accuracy for the same number of elements. However, creating the mesh is only one step in the simulation pipeline. The admittance method used here is only one approach in converting the mesh to a network of discrete conductances; it is straightforward and effective, but it is not necessarily the best pairing with this type of mesh. The next chapter will compare alternatives to this discretization method to see if its speed, accuracy, or both, can be exceeded. 54 Chapter 4 Electrode Representation and Activation Thresholds With confidence in the performance and scalability of this simulation pipeline, it can now be ap- plied to practical studies. The next chapter will move from the static cases covered thus far to fully time-varying sources; this will form the groundwork for accelerating a real-world application. First, however, this chapter will expand on the practical effects of electrodes’ meshed resolution on the neurons’ intracellular behavior. 4.1 Motivation Given its simplicity, approximating electrodes as point sources is still common in studies of intra- cellular behavior [1, 54, 27]. As distance from the source increases, and thus the geometric features of the electrode become proportionally smaller, the quality of this approximation should likewise increase—but its relevance to the neural tissue diminishes due to the decay in the field. For studies based on real, physical electrodes, accuracy will depend on a meshed representation of the actual electrode geometry. As is generally the case with these methods, it is difficult to establish a priori what level of meshing detail is necessary to attain a given degree of accuracy, beyond general heuristics like scaling elements proportionally to their distance from a source. With the wide range of resolutions possible in the octree meshing approach, a source’s rep- resentation may fall into one of two regimes: a single mesh node, or a collection of nodes that approximates the actual electrode geometry. Within the former regime, error within the adjacent 55 elements (compared against the analytical point source) can be several times the source’s magni- tude, especially when element size is close to the source radius. It is unclear whether the overshoot observed in the initial tests is likely to recur in a realistic simulation, thus requiring explicit safe- guards, or if it was simply an artifact of having a source at the domain’s origin. More generally, it would be useful to have a pre-simulation estimate of the expected detail required to accurately represent a stimulating electrode’s effect on the arrangement of neurons, if only to plan the computational resources it will require. If it is known in advance that a point- source approximation is sufficient, the added burden of a meshed simulation might be avoided; if a subset of neurons will require the explicit meshed representation, then at least the others can be safely omitted from evaluation of the mesh adaptation. 4.2 Initial study: toy neuron 4.2.1 Methods For initial proof of concept studies, a compartmental model with minimal geometric complexity is sufficient as long as its membrane dynamics reasonably capture the dynamics of a real neuron. Since the NEURON software package [22, 55] is thedefacto standard compartmental model sim- ulator, it will likewise be used here for all intracellular dynamics; specifically, the ball-and-stick neurons used in its tutorials[56] provide an appropriate toy model for preliminary investigations. Key parameters of this model are given in table 4.1. Due to the wide variance of reported conduc- tivities in the literature, an ensemble value ofσ = 0.3841S/m for whole-brain conductivity will be used for the extracellular space to provide a plausible scenario [42]. The stimulating electrodes are based on an insulated microwire pair from prior hippocampal studies [57, 26], whose conductive surfaces take the form of 48µm diameter disks on a 96µm pitch. However, since these electrodes are thus separated by a distance of only two electrode radii, it is impossible to vary the size of local elements without using a fine mesh (with multiple nodes within the electrode boundary), or significantly distorting the effective location of the electrodes 56 Table 4.1: Toy neuron parameters Region Quantity Value Unit Soma Diameter 12.6157 µm Dendrite Diameter 1 Length 200 Segments 5 All ρ axial 100 Ω· cm C m 1 µF cm 2 Soma ¯ g Na 0.12 S cm 2 ¯ g K 0.036 g l 0.0003 E l -54.3 mV Dendrite g p 0.001 µF cm 2 E p -65.0 mV by ”snapping” them to a distant node. Consequently, a single microwire will instead be used to deliver the stimulus, allowing a more granular sweep of the element geometry that is more likely to capture any overshoot effects that would significantly distort the apparent stimulus threshold. For the initial proof-of-concept tests, the stimulus takes the form of square biphasic current pulses at 1ms per phase, in line with typical stimulus protocols [57, 1, 58]. Since the initial elec- trode current i src (t 0 ) results in an extracellular potential v n (t 0 ) at neural compartment n, a single simulation on the mesh allows determination of an effective resistance R n = v n (t 0 ) i src (t 0 ) (4.1) that will be equal for all time points and possible current amplitudes, so the time course of the extracellular constraint is given by v n (t)=i src (t)R n (4.2) Note that the neural compartments do not necessarily need to correspond to mesh nodes, as their extracellular voltage can be interpolated from the nodal voltages of their containing element [24, 26]. The simulation mesh can thus be far sparser during the stimulus phase, provided it adequately models the high-gradient regions of the domain. 57 The change in this extracellular potential is either sufficient to generate an action potential or not; by starting with sufficiently strong stimulus, the minimum amplitude to elicit the spike can be found by binary search on the amplitude, re-performing the intracellular simulation at each tested amplitude. Because the overshoot phenomenon is most pronounced in the space between the source node and its immediate neighbors, more distant neurons should not be affected; thus, neurons placed up to 150 µm from the electrode center should be sufficient to delineate the region where neurons are susceptible to the effect. Figure 4.1: Examples of extracellular potential at peak of stimulus, with electrode placed 25µm (left) and 150µm (right) from neuron. The model neuron consists of a spherical soma with Hodgkin-Huxley dynamics and a passive dendrite. The cellular compartments are independent of the mesh, with their external potential determined by interpolation. 4.2.2 Results Surprisingly, variance in threshold as a function of resolution was most pronounced not in the nearby neurons, but the distant ones. Figure 4.2 shows relatively little change of threshold in the overshoot region compared to finer meshes, but the fine meshes (with multiple mesh nodes representing the source) exhibit large spikes in threshold only for the more distant neurons. Figure 4.3 shows these simulated thresholds normalized by the threshold of the equivalent analytic point-source stimulus. Neurons at all distances tend toward a ratio between 2 and 3, 58 Figure 4.2: Variation in activation threshold as a function of local mesh resolution. Significant variance occurs for more distant neurons in meshes finer than the overshoot region, where conver- gence would be expected. 59 albeit with significant variation; the more distant neurons in particular swing by nearly an order of magnitude. Figure 4.3: Ratio of activation thresholds for meshed disk electrode vs. analytic point-source. Though noise is substantial at all resolutions, neurons at all distances converge to a ratio of 2.5 for a fine mesh. A possible explanation for these inconsistencies is given in figure 4.4, where the gradient be- tween soma and dendrite varies more gradually in the near-field neuron than in the distant one. Since the gradient within an admittance-method element is constant, but no continuity is enforced between adjacent elements, there can be a significant difference in the voltage between two points depending on whether they lie within the same element or not. As the mesh resolution is swept, a distant neuron will ”pass through” more elements of differing sizes than one close to the source, and thus will experience more variance in external potential between simulations. 60 Figure 4.4: Extracellular potential of each compartment as mesh resolution varies, for neuron 50 (left) and 150µm (right) from electrode. Variance of gradient between soma (≈− 100µm) and remainder of the neuron at different resolutions is likely responsible for the wide variance in thresholds of the distant neuron. 4.3 Realistic neuron model The two-compartment toy model is clearly too simplistic to provide meaningful insights, and so a more realistic substitute is needed. Its geometry is too compact, leading to undue sensitivity to small changes in the external potential gradient; lacking an axon, it cannot be reasonably extended to the dimensions of a realistic neuron. Instead, a generic but realistic neuron model is called for. Figure 4.5 illustrates several common models from the literature. Though its protocol involved point electrodes at much further distances than will be studied here, it serves to illustrate the significant spread in excitability among models that intend to capture similar biology. Of these, the MRG model will be used due to its easy avail- ability and frequent adaptation in the literature [54, 1]. Because near-field effects are of interest in this study, a much smaller neuron of only 9 nodes of Ranvier is sufficient to extend well past the range of meaningful stimulus and avoid boundary errors. As in the previous study, the relative positions of electrode and neurons can be seen in figure 4.6, though the simulation domain extends much farther (and the axon further still) beyond the area shown. The activation thresholds are shown in figure 4.7, with much clearer trends than the results of the toy model. Again, the anticipated overshoot region does not exhibit any notable distortion 61 Figure 4.5: Comparison of strength-duration curves for common models of myelinated neurons, reproduced from [1]. 62 in the activation threshold; however, the variations in the fine meshes are much more subdued at all distances, and a fairly small percentage of their total. Compared to the equivalent analytic point-sources, thresholds are 2-5x higher, with closer neurons behaving more like the point-source solution. Figure 4.6: Detail of MRG neuron model 25µm (left) and 150µm (right) away from the electrode. 4.4 Discussion and conclusion Though the potential field in the immediate vicinity of a current source may be severely overesti- mated at certain resolutions, this phenomenon seems highly unlikely to affect any practical simula- tions. Because it occurs only at a narrow range of element sizes just larger than the electrode itself, it already requires an undesirable combination of relatively high element count that still reduces the source to a single node of the mesh. Furthermore, any other nearby sources (as in a bipolar pair) will either constrain the possible element sizes to multiples of their radii, or introduce other errors due to moving the relative location of the electrodes. Such a combination of scenarios could only arise rarely as an artifact of an unsupervised adaptation procedure, as any manually set-up simulation would involve much higher resolution to capture the electrode geometry. Even in this worst-case confluence of events, all the other sources of uncertainty contribute more to the effective threshold than this overestimation. Because it is such a local phenomenon, 63 10 0 10 1 ℓ 0 /r elec 50 μA 100 μA 150 μA 200 μA 250 μA 300 μA 350 μA 400 μA Threshold Overshoot region Activation threshold for 1.0 ms biphasic pulse 0μμ 25μμ 50μμ 75μμ 100μμ 125μμ 150μμ (a) 10 0 10 1 ℓ 0 /r elec 10 −1 10 0 10 1 10 2 i sim /i analytic Threshold ratios 0μm 25μm 50μm 75μm 100μm 125μm 150μm (b) Figure 4.7: Activation threshold of MRG neurons by mesh resolution and distance from electrode. (a) No difference in threshold between overshoot region and other mesh sizes. (b) Thresholds due to the modeled disk electrode converge to 2-3x the threshold of the analytic point source, except in the immediate vicinity of the electrodes. 64 most of the neuron remains unaffected even if it passes directly underneath the electrode, with negligible change to the net threshold. One unexpected finding is that the disk electrodes yield thresholds closer to those of a point source the closer they are to the neuron, whereas intuition suggests the field should become in- creasingly point-like at a distance from the electrode. While this is not a concern for meshed simulations, it does suggest that a correction factor may be required even for distant analytical sources if they are meant to represent a real-world geometry. 65 Chapter 5 Dynamic Mesh Adaptation While the earlier investigations kept to static conditions at a single point in time, any useful simu- lation will vary in time as well. In the literature, many simulations have used a single mesh for the entirety of the simulated time course, though the variation in activity patterns means it is far from optimal at any given point. If adaptation is itself computationally cumbersome, this mesh reuse may still be the most efficient compromise; there is little sense in adding an operation that takes longer than the time it saves in later steps, except perhaps when memory is the limiting factor and accuracy is paramount. As implemented, adapting the mesh occupies a significant fraction of the total simulation time—too long to ignore the computational requirements, yet short enough its cost is likely justi- fied by later savings. Where this break-even point lies will vary with the specifics of a simulation: how many sources are there, how widely their amplitude ranges, and how diffuse their spatial arrangement. To give a rough indication of how efficiency scales with these conditions, simulations on toy neural models are performed both on static meshes and dynamically adapted ones. A single source provides the first scenario tested, and a handful of generic neurons the second. The trend toward greater efficiency as complexity increases suggests this adaptive meshing approach will be a useful tool in real-world applications. 66 5.1 Single compartment The simplest use case for dynamic mesh adaptation is that of the single, stationary source whose amplitude varies over time. In this case, times of higher amplitude are expected to warrant higher mesh resolution, so the source’s absolute value is used to scale the effective resolution of the mesh. 5.1.1 Methods Again, the minimal toy neuron is used to generate the source waveform; the current through the soma’s membrane resulting from the artificial dendritic stimulus is used to drive a spherical source at the origin, the same radius as the soma (around 6 µm) It is undesirable to reduce mesh resolution too far even when amplitudes are particularly low, as this will still contribute noticeable error at these time points. Instead, it is convenient to define an approximate recursion depth at timet according to N ′ =n 0 +(n 1 − n 0 ) |i src (t)| max|i src | (5.1) such that element splitting is capped at⌊N ′ ⌋ and the static size field of equation 3.1 becomes ℓ 0 =∥r∥2 − kN ′ (5.2) Whereas earlier sections dealt strictly with increasing mesh resolution, a real simulation will experience ebb and flow in activity levels, not monotonic increase. To truly adapt, the mesh must thus have a means for pruning unnecessary detail as well as adding it when needed. One approach is to simply use the same size field that drives refinement; instead of splitting whenever l>ℓ 0 , an octant can safely purge its children if l<ℓ 0 for itself and all the children. In this way, the octant is now a terminal leaf in the octree, and replaces its children in the mesh. Even using the same function for both refinement and pruning, this method introduces a degree of hysteresis, and thus 67 avoids unnecessary cycling between the two. As implemented, a different size field can instead be used for pruning if the user desires, defaulting to the refinement field otherwise. 5.1.2 Results The resulting potential field of the basic adapted mesh, at a representative set of frames, is shown in figure 5.1. Figure 5.2 gives a detailed comparison of the accuracy and simulation time per timestep between the adapted mesh, and static octree meshes equivalent to maximum and minimum adapted resolution. A clearer picture emerges in figure 5.3 showing the per-simulation totals of error and execution time. 68 Figure 5.1: Images of single-source simulation with dynamically adapted mesh 69 In this single-source scenario, adaptation is not always the hands-down best performer; rather, it provides greater control in the tradeoff between accuracy and speed. It is essentially a midpoint between the minimal- and maximal-resolution static meshes: around half the error in twice the time compared to the minimal mesh, and a fraction the time of the max-resolution mesh at far lower accuracy. Allowing even deeper recursion resulted in similar total execution times as the fixed maximal-resolution mesh, and failed to match its accuracy though it significantly improved on the shallower adaptation. Figure 5.2: Accuracy and computational speed of fixed and dynamically adapted meshes. 70 Figure 5.3: Total error and computational time for fixed and dynamic meshes 71 Though a needed proof-of-concept, this experimental setup is not a reasonable choice for dy- namic adaptation. With a single source, it is far faster to simulate a high-resolution mesh once to find v(t 0 ), and calculate the subsequent timesteps by v(t)= v(t 0 ) i src (t) i src (t 0 ) (5.3) Though that method only applies to a single source, it can still be an economical approach for a small number of sources if their individual fields are added by superposition. 5.2 Toy multi-neuron model 5.2.1 Methods Multiple toy neurons in a ring pattern are now used for current sources in the extracellular simula- tion, with an artificial synapse from the soma of each to the dendrite of its neighbor. The frames of figure 5.4 indicate this layout as well as the resulting field. The artificial stimulus described in the previous section triggers the initial neuron’s action potential, after which all electrical activity is due to the neurons’ internal biophysics. 40ms of transmembrane current is recorded from every neural compartment and injected at the corresponding mesh location. Since there are now multiple sources in the mesh, calculation of ℓ 0 is a function of distance to each source. Though each could conceivably affect the desired resolution, most will of course be too far away to matter. Thus, at each step of testing the mesh refinement, only the calculations that would split a parent element is passed on to its children, discarding those that will also be irrelevant within the child octants. Likewise, the desired depth is scaled according to the maximum of all sources, but calculated per-source (and filtered to children) as as described above. 72 5.2.2 Results As before, figure 5.4 contains slice images of the resulting fields, figure 5.5 the error and computa- tion time per time step, with the per-simulation totals shown in figure 5.6. Unlike the single-source case, there is a dramatic separation in simulation times between the adaptive cases, which are nearly as fast as the low-resolution mesh, and the fixed high-resolution mesh, which is has error equivalent to the much faster adaptive meshes. Unfortunately, increasing the recursion depth only produced a moderate reduction in total error after doubling execution time. In other words, it seems the benefits of this adaptation process are better described in terms of the reduction of execution time for equivalent accuracy, rather than of increased accuracy per resources; at the tested maximum of 12 generations, the smallest element is already around 300 nm per side, and well past the point of diminishing returns. 5.3 Realistic model geometry The improvements seen in the toy multi-neuron model suggest benefits when there are only a handful of independently-active sources, as might occur in a multi-electrode stimulating array. It is expected to be even more beneficial in estimating local field potentials in a large-scale neural model, due to the sheer number of independent sources involved. However, implementing such a model will require significant resources for the intracellular domain, independent of the extra- cellular methods employed here. Absent a ready-made dataset of neural compartments’ locations and transmembrane currents, a more tractable demonstration of the dynamic adaptation process is in order; hence, a stimulating multielectrode array within realistic tissue geometries provides a logical final demonstration of dynamic adaptation. 5.3.1 Methods The simulation is based on prior in-vivo studies in the human hippocampus [58]. Figure 5.7 gives an overview of the experimental setup: a commercial DBS electrode array is implanted in the 73 Figure 5.4: Images of simulation with dynamically adapted mesh 74 Figure 5.5: Accuracy and computational speed of fixed and dynamically adapted meshes. 75 Figure 5.6: Total error and computational time for fixed and dynamic meshes hippocampus, with the microelectrodes, represented by the gold dots between the larger band macroelectrodes, present in both CA1 and CA3. Tissue regions within the hippocampus can have substantially different conductivities [59, 57]. Though 3-dimensional anatomical datasets are readily available for human anatomy, most do not have sufficient resolution to reliably distinguish these different substructures. Instead, a high- resolution 2d atlas forms the basis of the hippocampal model [60]. After combining regions with the same conductivity (per the schema in [57]), their net boundary is extruded to define its 3d boundary. Following [58], a single stimulus consists of a biphasic square pulse with 150µA amplitude and 1ms per phase delivered to one of the microelectrodes. The overall pattern of stimulation is governed by pseudorandom Poisson processes averaging 10Hz per channel, binned at 50ms intervals to ensure a minimal inter-pulse duration. In the present study, all 24 microelectrode channels are used. 76 (a) (b) Figure 5.7: Schematic of DBS electrode in hippocampus (a) and slice of minimal resolution mesh (b). 77 Each microelectrode was first tested independently to determine the mesh resolution required for convergence of the solution; all were found to change minimally after 16 subdivisions. Con- versely, a minimum of 6 subdivisions barely resolves the geometric features of the hippocampus, as seen in figure b, and serves as a lower bound on resolution for the adaptive process. To keep the total simulation times more tractable, the density factor was reduced from the 0.2 used in prior studies to 0.1. The geometric boundaries are not considered while constructing the mesh, only the location (and in the dynamic case, amplitude) of the sources. Instead, the region to which each element belongs is set by which original surface mesh contains its center. For the sake of development expedience and to guarantee every element expressed the correct conductivity, for each timestep where the mesh topology changes the elements are re-assigned region membership based on these high-resolution surface meshes. In this way, the children of an element that crosses a region bound- ary can receive the conductivity appropriate to which side they lie on, regardless of the parent’s conductivity. Unlike the previous studies, the heterogeneous conductivities make it impossible to calculate an equivalent analytical solution, and therefore an error estimate. Instead, proof of accuracy must be established by the difference between the solution on a static mesh and the solution on a dynamically-adapted one, not between either and an analytical function. Since individual stim- uli do not change appreciably for meshes beyond 16 divisions, it seems safe to treat the static, maximally-refined mesh as a proxy for the ground truth in this fashion, and evaluate the dynamic adaptation based on its time savings versus deviation in the solution. 5.3.2 Results Figure 5.8 gives a representative sample of the simulated field under dynamic adaptation. The electrodes intersected by the image plane are easily distinguished from those out of the plane by their stronger field; the presence of the macroelectrodes’ equipotential surface can also be seen in the interior of the probe body. 78 The difference in solutions between the static and dynamically adapted meshes is at most a few percent, as shown by figure a. Nonetheless, the dynamic simulation is around an order of magnitude faster, as seen in figure b, despite needing to repeat several stages of the meshing process at every timestep. 5.4 Discussion and conclusions In simulations of both local field potential and external electrical stimulus, dynamic mesh adap- tation can reduce simulation times by an order of magnitude while maintaining accuracy of the results. Because the spatial distribution of the field varies so much over the course of the simula- tion, detail in a given region of the domain will only contribute to solution accuracy in a fraction of the timesteps, while only inflating the computational burden in the rest of the simulation. Though the final step of solving the system GV =I can be efficiently handled by mature, well- optimized iterative routines, it is still a computationally intensive process that must be repeated at every timestep, regardless of what else may occur in the simulation. Consequently, the most tractable way of accelerating it (without adopting more powerful hardware) is to reduce the size of the system, i.e. minimizing the number of elements in the mesh. In contrast, other time-consuming steps of the current implementation (like the node sorting that dominates in figure b) are mostly custom Python routines that could likely be reimplemented as vastly faster optimized and compiled code. If so, the adaptation process would prove even more advantageous, as the final solving step becomes a greater fraction of the solution time. Similarly, the hippocampal stimulation test in some ways understates the advantages of the adaptation process. Since its stimuli are essentially binary (either full amplitude or quiescent at a given timestep), the limitation of local pruning to a single level per timestep means detail will per- sist several times longer than needed to adequately model the pulse. A model that incorporates the capacitance of the electrode-tissue interface, or otherwise has a finer time-domain representation, would result in a closer match between the local potential amplitude and the actual resolution. 79 Figure 5.8: Potential field in the hippocampus under pseudorandom stimulation. 80 (a) (b) Figure 5.9: Summary comparison of fixed and dynamically-adapted meshes in the hippocampal simulation. (a) Though much faster than the fixed, high-resolution mesh, the field simulated by the dynamic mesh differs by only a few percent. (b) As in the LFP estimation, net simulation time is an order of magnitude faster under dynamic adaptation. Despite the inefficiency of re-assigning conductivity at every timestep, it is responsible for a small portion of the total time spent in the dynamic adaptation. 81 Chapter 6 Conclusions To summarize, this dissertation seeks to make four main points: • Mathematical basis of mesh adaptation: when treating the extracellular space as a resistor network, it is impossible to locally adapt the topology while maintaining electrical equiva- lence. Rather, a network approximation should be applied at each geometric region. • Octrees guide efficient meshing: by recursively subdividing, a mesh can be quickly con- structed, with a spatial resolution tailored to local requirements. Thus, computational re- sources are efficiently distributed, outperforming the naive approach of a uniform grid; con- versely, the simplest method of approximating local conductance is also the most accurate when applied to these sort of meshes. • Point current sources can cause overshoot: because the potential of a current source ap- proaches infinity as its size approaches a point, there can be significant error introduced at mesh resolutions near the source radius. However, because the effect is so localized, it has minimal practical impact on activation thresholds. • Dynamic mesh adaptation further speeds simulation: if there are more than a handful of independent sources in a simulation, distributed throughout the domain, the overhead of remeshing when the locus of activity changes is far outweighed by the reduction in the size of the system of equations. 82 Though the preliminary exploration of dynamic adaptation in this work shows its potential benefits in neural simulations, it is still a preliminary effort. The combination of the admittance method with an underlying octree structure forms a solid basis for a robust simulation pipeline, but there are still unexplored improvements both in optional stages of the simulation process, and in the concrete algorithmic implementations of those already demonstrated. Elements need not be split evenly, nor does the a priori estimate of mesh resolution need to be its sole determinant. Adjusting the proportion of an element split allows placement of a node at any location within the parent element without an increase in the size of the overall system; in broader finite-element literature, this approach is known as r-refinement, frequently employed in fluid dynamics simulations [61, 62]. In the present context, this node-shifting could allow better representation of geometric boundaries prior to solving for the potential field, as well as the traditional a posteriori application in re-simulating a timestep for improved accuracy. Of course, any error metric used to drive this refinement could as easily be applied to the splitting-based h- refinement demonstrated in this work, and further experimentation is needed to indicate the relative merits of each approach. Though equivalent-circuit transforms are ruled out as impossible, and mesh-dual formulations as computationally inefficient, there may yet be value in a sort of transform applied to the electrical network of the simulation. In theory, the changes in the mesh due to refinement could be isolated, and defined in terms of the original mesh quantities; with these substitutions cast as matrix opera- tions, as in some adaptation approaches [63], there may be further time savings through avoidance of re-calculating quantities or simply taking advantage of extant optimizations for matrix routines. Similarly, there are surely optimizations possible simply from restructuring the methods em- ployed here. As implemented, the meshing processes are defined within the context of the element, and so every split could conceivably generate 8 independent threads until all available processors are in use. Due to its added complexity and potential for error, however, this multithreading was not yet implemented. An alternative approach to parallelism would be to treat all elements of the same 83 depth together, in a single-instruction multi-data (SIMD) operation. This paradigm might ironi- cally be more flexible, due to its prevalence in GPU programming and other numerically intense (but essentially repetitive) operations. Finally, all efforts at improving simulations of electrode-neuron interactions are useless without robust in-vitro and in-vivo experimentation to validate their predictions and justify their assump- tions. Otherwise, as model complexity grows to take advantage of computational advances, it can create the illusion of insight without basis in reality. Conversely, proliferation of powerful yet easy-to-use simulation tools can be a tremendous boon to the experimentalist, allowing rapid pro- totyping of hypotheses before and throughout the real-world data collection. Public release of the xcell library developed for this dissertation, along with its associated documentation, hopes to be a small contribution toward this goal. 84 References [1] J. P. Reilly, “Survey of numerical electrostimulation models,” Phys. Med. Biol., vol. 61, pp. 4346–4363, May 2016. Publisher: IOP Publishing. [2] T. W. Berger, D. Song, R. H. M. Chan, V . Z. Marmarelis, J. LaCoss, J. Wills, R. E. Hamp- son, S. A. Deadwyler, and J. J. Granacki, “A Hippocampal Cognitive Prosthesis: Multi-Input, Multi-Output Nonlinear Modeling and VLSI Implementation,” IEEE Transactions on Neu- ral Systems and Rehabilitation Engineering, vol. 20, pp. 198–211, Mar. 2012. Conference Name: IEEE Transactions on Neural Systems and Rehabilitation Engineering. [3] J. D. Weiland and M. S. Humayun, “Retinal Prosthesis,” IEEE Transactions on Biomedical Engineering, vol. 61, pp. 1412–1424, May 2014. Conference Name: IEEE Transactions on Biomedical Engineering. [4] S. Elyahoodayan, C. Larson, A. M. Cobo, E. Meng, and D. Song, “Acute in vivo testing of a polymer cuff electrode with integrated microfluidic channels for stimulation, recording, and drug delivery on rat sciatic nerve,” Journal of Neuroscience Methods, vol. 336, p. 108634, Apr. 2020. [5] A. L. Hodgkin and A. F. Huxley, “A quantitative description of mem- brane current and its application to conduction and excitation in nerve,” The Journal of Physiology, vol. 117, no. 4, pp. 500–544, 1952. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1113/jphysiol.1952.sp004764. [6] H. Lind´ en, E. Hagen, S. Leski, E. S. Norheim, K. H. Pettersen, and G. T. Einevoll, “LFPy: a tool for biophysical simulation of extracellular potentials generated by detailed model neu- rons,”Front.Neuroinform., vol. 7, 2014. Publisher: Frontiers. [7] C. H. Lubba, Y . Le Guen, S. Jarvis, N. S. Jones, S. C. Cork, A. Eftekhar, and S. R. Schultz, “PyPNS: Multiscale Simulation of a Peripheral Nerve in Python,” Neuroinform, vol. 17, pp. 63–81, Jan. 2019. [8] M. Capllonch-Juan and F. Sepulveda, “Evaluation of a Resistor Network for Solving Electri- cal Problems on Ohmic Media,” in201911thComputerScienceandElectronicEngineering (CEEC), pp. 35–40, Sept. 2019. [9] S. Appukuttan, K. L. Brain, and R. Manchanda, “Modeling extracellular fields for a three- dimensional network of cells using NEURON,” Journal of Neuroscience Methods, vol. 290, pp. 27–38, Oct. 2017. 85 [10] A. P. Buccino, M. Kuchta, K. H. Jæger, T. V . Ness, P. Berthet, K.-A. Mardal, G. Cauwen- berghs, and A. Tveito, “How does the presence of neural probes affect extracellular poten- tials?,”J.NeuralEng., vol. 16, p. 026030, Feb. 2019. Publisher: IOP Publishing. [11] S. Joucla, A. Gli` ere, and B. Yvert, “Current approaches to model extracellular electrical neural microstimulation,”Front.Comput.Neurosci., vol. 8, 2014. [12] A. Agudelo-Toro and A. Neef, “Computationally efficient simulation of electrical activity at cell membranes interacting with self-generated and externally imposed electric fields,” J. NeuralEng., vol. 10, p. 026019, Mar. 2013. Publisher: IOP Publishing. [13] A. P. Buccino, M. Kuchta, J. Schreiner, and K.-A. Mardal, “Improving Neural Simulations with the EMI Model,” in Modeling Excitable Tissue (A. Tveito, K.-A. Mardal, and M. E. Rognes, eds.), vol. 7, pp. 87–98, Cham: Springer International Publishing, 2021. Series Title: Simula SpringerBriefs on Computing. [14] Y . Hu, Q. Zhou, X. Gao, A. Jacobson, D. Zorin, and D. Panozzo, “Tetrahedral Meshing in the Wild,” ACM Trans. Graph., vol. 37, July 2018. Place: New York, NY , USA Publisher: Association for Computing Machinery. [15] X. Gao, W. Jakob, M. Tarini, and D. Panozzo, “Robust hex-dominant mesh generation using field-guided polyhedral agglomeration,” ACMTrans.Graph., vol. 36, pp. 114:1–114:13, July 2017. [16] X. Gao, H. Shen, and D. Panozzo, “Feature Preserving Octree-Based Hexahedral Meshing,” Computer Graphics Forum, vol. 38, no. 5, pp. 135–149, 2019. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.13795. [17] T. Schneider, Y . Hu, X. Gao, J. Dumas, D. Zorin, and D. Panozzo, “A Large-Scale Com- parison of Tetrahedral and Hexahedral Elements for Solving Elliptic PDEs with the Finite Element Method,”ACMTrans.Graph., vol. 41, pp. 23:1–23:14, Mar. 2022. [18] C. J. Cela, A Multiresolution Admittance Method for Large-Scale Bioelectromagnetic Inter- actions. PhD Thesis, North Carolina State University, 2010. ISBN: 978-1-124-76976-9 Publication Title: ProQuest Dissertations and Theses. [19] R. Duffin, “Distributed and Lumped Networks,” IndianaUniv.Math.J., vol. 8, no. 5, pp. 793– 826, 1959. [20] M. Komarudin, “Resistor Network Analogy of Non-obtuse Finite Element Model for Elec- trical Impedance Tomography,”Electrician, vol. 1, pp. 1–10, Sept. 2007. Number: 1. [21] A. Al-Humaidi, Resistor Networks And Finite Element Models. PhD Thesis, University of Manchester, 2011. ISBN: 9781073271948 Publication Title: PQDT - UK & Ireland. [22] M. Hines, A. Davison, and E. Muller, “NEURON and Python,” Frontiers in Neuroinformat- ics, vol. 3, 2009. 86 [23] P. J. Hendrickson, G. J. Yu, D. Song, and T. W. Berger, “A Million-Plus Neuron Model of the Hippocampal Dentate Gyrus: Critical Role for Topography in Determining Spatiotemporal Network Dynamics,” IEEE Transactions on Biomedical Engineering, vol. 63, pp. 199–209, Jan. 2016. Conference Name: IEEE Transactions on Biomedical Engineering. [24] P. J. Hendrickson, C. Bingham, D. Song, and T. W. Berger, “A bi-directional communication paradigm between parallel NEURON and an external non-neuron process,” in 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pp. 1413–1416, Aug. 2016. ISSN: 1558-4615. [25] A. Fellner, A. Heshmat, P. Werginz, and F. Rattay, “A finite element method framework to model extracellular neural stimulation,” J. Neural Eng., vol. 19, p. 022001, Apr. 2022. Publisher: IOP Publishing. [26] C. S. Bingham, J. Paknahad, C. B. C. Girard, K. Loizos, J.-M. C. Bouteiller, D. Song, G. Lazzi, and T. W. Berger, “Admittance Method for Estimating Local Field Potentials Gen- erated in a Multi-Scale Neuron Model of the Hippocampus,” Frontiers in Computational Neuroscience, vol. 14, 2020. [27] C. A. Bossetti, M. J. Birdno, and W. M. Grill, “Analysis of the quasi-static approximation for calculating potentials generated by neural stimulation,”J.NeuralEng., vol. 5, pp. 44–53, Mar. 2008. [28] K. Truemper, “On the delta-wye reduction for planar graphs,” Jour- nal of Graph Theory, vol. 13, no. 2, pp. 141–148, 1989. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/jgt.3190130202. [29] L. Versfeld, “Remarks on star-mesh transformation of electrical networks,” Electronics Let- ters, vol. 6, no. 19, p. 597, 1970. [30] V . V . B. Rao and V . K. Aatre, “Mesh-star transformation,”ElectronicsLetters, vol. 10, pp. 73– 74, Mar. 1974. [31] C. J. Cela, R. C. Lee, and G. Lazzi, “Modeling Cellular Lysis in Skeletal Muscle Due to Electric Shock,” IEEE Transactions on Biomedical Engineering, vol. 58, pp. 1286–1293, May 2011. Conference Name: IEEE Transactions on Biomedical Engineering. [32] C. Geuzaine and J.-F. Remacle, “Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities,” International Journal for Numer- ical Methods in Engineering, vol. 79, no. 11, pp. 1309–1331, 2009. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.2579. [33] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Po- lat, Y . Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,”NatureMethods, vol. 17, pp. 261–272, 2020. 87 [34] Z. Xiong, S. Feng, R. Kautz, S. Chandra, N. Altunyurt, and J. Chen, “Multi-GPU Accelerated Admittance Method for High-Resolution Human Exposure Evaluation,” IEEE Transactions on Biomedical Engineering, vol. 62, pp. 2920–2930, Dec. 2015. Conference Name: IEEE Transactions on Biomedical Engineering. [35] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del R´ ıo, M. Wiebe, P. Peterson, P. G´ erard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” Nature, vol. 585, pp. 357–362, Sept. 2020. Number: 7825 Publisher: Nature Publishing Group. [36] E. van der Velden, “CMasher: Scientific colormaps for making accessible, informative and ’cmashing’ plots,” TheJournalofOpenSourceSoftware, vol. 5, p. 2004, Feb. 2020. eprint: 2003.01069. [37] W. McKinney, “Data Structures for Statistical Computing in Python,” in Proceedings of the 9thPythoninScienceConference (S. v. d. Walt and J. Millman, eds.), pp. 56 – 61, 2010. [38] J. D. Hunter, “Matplotlib: A 2D graphics environment,” Computing in Science & Engineer- ing, vol. 9, no. 3, pp. 90–95, 2007. Publisher: IEEE COMPUTER SOC. [39] O. Awile, P. Kumbhar, N. Cornu, S. Dura-Bernal, J. G. King, O. Lupton, I. Magkanaris, R. A. McDougal, A. J. H. Newton, F. Pereira, A. S˘ avulescu, N. T. Carnevale, W. W. Lytton, M. L. Hines, and F. Sch¨ urmann, “Modernizing the NEURON Simulator for Sustainability, Portability, and Performance,”FrontiersinNeuroinformatics, vol. 16, 2022. [40] S. K. Lam, A. Pitrou, and S. Seibert, “Numba: a LLVM-based Python JIT compiler,” in Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC - LLVM ’15, (Austin, Texas), pp. 1–6, ACM Press, 2015. [41] C. B. Sullivan and A. A. Kaszynski, “PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK),”JournalofOpenSourceSoftware, vol. 4, p. 1450, May 2019. [42] H. McCann, G. Pisano, and L. Beltrachini, “Variation in Reported Human Head Tissue Elec- trical Conductivity Values,”BrainTopogr, vol. 32, pp. 825–858, Sept. 2019. [43] K. Muramatsu, T. Nakata, N. Takahashi, and K. Fujiwara, “Investigation of effectiveness of 3-D nonconforming mesh,” IEEE Transactions on Magnetics, vol. 27, pp. 5211–5213, Nov. 1991. Conference Name: IEEE Transactions on Magnetics. [44] S. Noguchi, T. Naoe, H. Igarashi, S. Matsutomo, and V . Cingoski, “A new adaptive meshing method using non-conforming finite element method,” in 2016IEEEConferenceonElectro- magneticFieldComputation(CEFC), pp. 1–1, Nov. 2016. [45] H. G. Kollech and J. P. Vande Geest, “Conforming versus nonconforming meshes in the finite element modeling of the human lamina cribrosa,” InvestigativeOphthalmology&Visual Science, vol. 62, p. 1643, June 2021. 88 [46] A. Gilbert, K. Loizos, A. K. RamRakhyani, P. Hendrickson, G. Lazzi, and T. W. Berger, “A 3-D admittance-level computational model of a rat hippocampus for improving prosthetic design,” in201537thAnnualInternationalConferenceoftheIEEEEngineeringinMedicine andBiologySociety(EMBC), pp. 2295–2298, Aug. 2015. ISSN: 1558-4615. [47] K. Loizos, G. Lazzi, J. S. Lauritzen, J. Anderson, B. W. Jones, and R. Marc, “A multi-scale computational model for the study of retinal prosthetic stimulation,” in201436thAnnualIn- ternationalConferenceoftheIEEEEngineeringinMedicineandBiologySociety, pp. 6100– 6103, Aug. 2014. ISSN: 1558-4615. [48] J. Paknahad, P. Kosta, J.-M. C. Bouteiller, M. S. Humayun, and G. Lazzi, “Mechanisms underlying activation of retinal bipolar cells through targeted electrical stimulation: a com- putational study,”J.NeuralEng., vol. 18, p. 066034, Dec. 2021. [49] B. Wang, A. Petrossians, and J. D. Weiland, “Reduction of Edge Effect on Disk Electrodes by Optimized Current Waveform,” IEEE Transactions on Biomedical Engineering, vol. 61, pp. 2254–2263, Aug. 2014. Conference Name: IEEE Transactions on Biomedical Engineer- ing. [50] B. Wang and J. D. Weiland, “Resistivity profiles of wild-type, rd1, and rd10 mouse retina,” in 2015 37th Annual International Conference of the IEEE Engineering in Medicine and BiologySociety(EMBC), pp. 1650–1653, Aug. 2015. ISSN: 1558-4615. [51] N. A. Pelot, B. J. Thio, and W. M. Grill, “Modeling Current Sources for Neural Stimulation in COMSOL,”Front.Comput.Neurosci., vol. 12, 2018. [52] J. E. B. Randles, “Kinetics of rapid electrode reactions,” Discuss. Faraday Soc., vol. 1, pp. 11–19, Jan. 1947. Publisher: The Royal Society of Chemistry. [53] Y . Hu, T. Schneider, B. Wang, D. Zorin, and D. Panozzo, “Fast tetrahedral meshing in the wild,”ACMTrans.Graph., vol. 39, pp. 117:117:1–117:117:18, July 2020. [54] C. C. McIntyre, A. G. Richardson, and W. M. Grill, “Modeling the Excitability of Mammalian Nerve Fibers: Influence of Afterpotentials on the Recovery Cycle,” JournalofNeurophysiol- ogy, vol. 87, pp. 995–1006, Feb. 2002. [55] N. T. Carnevale and M. L. Hines,TheNEURONbook. Cambridge University Press, 2006. [56] “Python tutorials — NEURON documentation.” [57] C. S. Bingham, K. Loizos, G. J. Yu, A. Gilbert, J.-M. C. Bouteiller, D. Song, G. Lazzi, and T. W. Berger, “Model-Based Analysis of Electrode Placement and Pulse Amplitude for Hip- pocampal Stimulation,” IEEE Transactions on Biomedical Engineering, vol. 65, pp. 2278– 2289, Oct. 2018. Conference Name: IEEE Transactions on Biomedical Engineering. [58] R. E. Hampson, D. Song, B. S. Robinson, D. Fetterhoff, A. S. Dakos, B. M. Roeder, X. She, R. T. Wicks, M. R. Witcher, D. E. Couture, A. W. Laxton, H. Munger-Clary, G. Popli, M. J. Sollman, C. T. Whitlow, V . Z. Marmarelis, T. W. Berger, and S. A. Deadwyler, “Developing a hippocampal neural prosthetic to facilitate human memory encoding and recall,” J. Neural Eng., vol. 15, p. 0360141v, Mar. 2018. Publisher: IOP Publishing. 89 [59] L. L´ opez-Aguado, J. M. Ibarz, and O. Herreras, “Activity-dependent changes of tissue resis- tivity in the CA1 region in vivo are layer-specific: modulation of evoked potentials,” Neuro- science, vol. 108, pp. 249–262, Dec. 2001. [60] S.-L. Ding, J. J. Royall, S. M. Sunkin, L. Ng, B. A. Facer, P. Lesnar, A. Guillozet-Bongaarts, B. McMurray, A. Szafer, T. A. Dolbeare, A. Stevens, L. Tirrell, T. Benner, S. Caldejon, R. A. Dalley, N. Dee, C. Lau, J. Nyhus, M. Reding, Z. L. Riley, D. Sandman, E. Shen, A. van der Kouwe, A. Varjabedian, M. Write, L. Zollei, C. Dang, J. A. Knowles, C. Koch, J. W. Phillips, N. Sestan, P. Wohnoutka, H. R. Zielke, J. G. Hohmann, A. R. Jones, A. Bernard, M. J. Hawrylycz, P. R. Hof, B. Fischl, and E. S. Lein, “Comprehensive cellular-resolution atlas of the adult human brain,”JournalofComparativeNeurology, vol. 524, no. 16, pp. 3127–3481, 2016. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/cne.24080. [61] D. S. McRae, “r-Refinement grid adaptation algorithms and issues,” Computer Methods in AppliedMechanicsandEngineering, vol. 189, pp. 1161–1182, Sept. 2000. [62] J. R. Grisham, N. Vijayakumar, G. Liao, B. H. Dennis, and F. K. Lu, “A Comparison Between Local h-Refinement and a Novel r-Refinement Method,” in 53rd AIAA Aerospace Sciences Meeting, (Kissimmee, Florida), American Institute of Aeronautics and Astronautics, Jan. 2015. [63] R. Anderson, J. Andrej, A. Barker, J. Bramwell, J.-S. Camier, J. Cerveny, V . Dobrev, Y . Du- douit, A. Fisher, T. Kolev, W. Pazner, M. Stowell, V . Tomov, J. Dahm, D. Medina, and S. Zampini, “MFEM: a modular finite element methods library,” Computers & Mathemat- icswithApplications, p. S0898122120302583, July 2020. arXiv: 1911.09220. 90 Appendices A Size field parameters To prevent infinite recursion nearest the source, the present method adds an explicit cap on recur- sion depthN. The two are combined into a single size field ℓ 0 =∥r∥2 − kN (A.1) where the parameter k roughly corresponds to the overall density of the mesh, as seen at the pe- riphery of the simulation domain; atk= 0, only the elements at the source center For a simulation domain 2x 0 along each axis, an element formed afterN subdivisions will have an edge length of l=x 0 2 1− N (A.2) and thus for the element with a vertex at the source, its center is ∥r∥= s 3 l 2 2 (A.3) = √ 3 2 l (A.4) 91 away from the source’s center. Splitting occurs ifl>ℓ 0 , so ∥r∥2 − kN <l (A.5) √ 3 2 l2 − kN <l (A.6) − kN< log 2 2 √ 3 (A.7) kN> log 2 2 √ 3 ≈ 0.2075 (A.8) At the other extreme, the furthest possible distance between source and element occurs for a source at one corner of the domain and the element at the opposite corner, where ∥r∥= q 3(2x 0 − l) 2 (A.9) = √ 3(2x 0 − l) (A.10) = √ 3(2x 0 − x 0 2 1− N ) (A.11) = 2 √ 3x 0 (1− 2 − N ) (A.12) and thus every element in the domain will be fully subdivided if ∥r∥2 − kN <l (A.13) 2 √ 3x 0 (1− 2 − N )2 − kN <x 0 2 1− N (A.14) √ 3(1− 2 − N )2 − kN < 2 − N (A.15) (1− 2 − N )2 N(1− k) < 1 √ 3 (A.16) 2 N(1− k) − 2 − k < 1 √ 3 (A.17) resulting in a uniform mesh resolution throughout. Substitutingk= 1, 92 2 N(1− k) − 2 − k < 1 √ 3 (A.18) − 1 2 < 1 √ 3 (A.19) and kN> log 2 2 √ 3 ≈ 0.2075 (A.20) N> 0.2... (A.21) which is valid for all values ofk andN, so the criteria are met for the entire mesh to be subdivided completely. At the other extreme, fork= 0 2 N(1− k) − 2 − k < 1 √ 3 (A.22) 2 N − 1< 1 √ 3 (A.23) (A.24) valid only forN= 0, i.e. no splitting of the mesh. Instead, from A.8 kN> log 2 2 √ 3 ≈ 0.2075 (A.25) N> 0.2... (A.26) 93 B Derivation of finite element method Formally, volume conductance is governed by ∇ 2 σV = f (B.1) where f is zero everywhere in the simulation domain except for at current sources/sinks (such that the nonzero value captures the current coming from ”outside” the domain). This equation is converted to the weak formulation by introducing another functionu so that Z Ω ∆ 2 σVudA= Z Ω fudA (B.2) where R Ω dA indicates integrating over the entire spatial domain (e.g., RRR dxdydz in 3d space). Applying Green’s identities yields − Z Ω ∇σV· ∇udA= Z Ω fudA (B.3) which sets us up for a final round of substitutions in terms of calculable quantities. The standard approach (known as Galerkin’s method) considers these quantities in terms of basis functions—usually, the sort referred to as ”hat” functions. Essentially, a hat function is 1 at a given node, zero at all other nodes, and smoothly varies between the given node and its immediate neighbors. Thus, each hat interpolates the value of the potential across a small region, and their sum covers the entire simulation domain. Essentially, each ”splits” the domain into small geometric sections, where relatively simple formulae can approximate the dynamics within that region. In the following discussions,h i will denote the hat function pertaining to nodei. 94 In this framework, we can call the voltage at node i v i , and the approximate solution of the overall potential field V ≈ ∑ N− 1 i=0 v i h i for all N nodes. We can further use these hat functions in place ofu in the earlier equation. With those substitutions, we arrive at − N− 1 ∑ i=0 v i Z Ω σ∇h i · ∇h j dA= Z Ω fh i dA (B.4) for each node. f will be zero, except in cases where current enters or leaves the simulated domain; likewise, the hat functions are zero except for inside elements containing the associated node. The hat functions depend only on mesh geometry, so all quantities can be directly calculated except for our unknownv i . We finally have a simple matrix equation Gv= i, with v a vector of the N v i unknowns, i the right-hand side of the equation, and element g i,j = Z Ω σ∇h i · ∇h j dA (B.5) 95

Abstract (if available)

Abstract While compartmental models have scaled to millions of neurons across thousands of cores, high resource consumption and limited parallelism has restricted simulation of the extracellular domain to comparatively low-fidelity approximations. Moreover, neural activity exhibits significant spatiotemporal variation, so reallocating detail at each timestep could speed simulation and improve accuracy.

This adaptation is considered as both an equivalent-circuit transform on the electrical network, and as an independent manipulation of mesh geometry. The former involves mutually incompatible constraints and is an unsuitable theoretical basis; the latter approach instead yields reasonable approximations.

The result is a simple, easily parallelized, means of discretizing the extracellular space according to activity. By recursively splitting until a target size is reached, the entire mesh maintains an octree structure amenable to localized alteration. With element size directly proportional to distance from the source, octree meshes achieve accuracy similar to a uniform grid of their smallest element, at a fraction of the computational burden.

Point sources reveal a potential pitfall as element sizes approach the source's radius, where potential in central elements significantly overshoots before additional mesh nodes enter the source and begin explicitly representing the boundary. Though naive simulation of injected current could thus yield inaccurate activation of nearby neural compartments, in practice the highly localized phenomenon does not affect overall excitability.

Finally, dynamic mesh adaptation yields accuracy comparable to a single, fixed mesh over the course of local field potential and multi-electrode stimulus simulations, but is an order of magnitude faster.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

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

Design and application of a C-shaped miniaturized coil for transcranial magnetic stimulation in rodents

PDF

Functional consequences of network architecture in rat hippocampus: a computational study

PDF

Multiscale spike-field network causality identification

PDF

Anatomically based human hand modeling and simulation

PDF

Homogenization procedures for the constitutive material modeling and analysis of aperiodic micro-structures

PDF

Electrical stimulation approaches to restoring sight and slowing down the progression of retinal blindness

PDF

Differential verification of deep neural networks

PDF

Development and evaluation of high‐order methods in computational fluid dynamics

PDF

A million-plus neuron model of the hippocampal dentate gyrus: role of topography, inhibitory interneurons, and excitatory associational circuitry in determining spatio-temporal dynamics of granul...

PDF

Neural spiketrain decoder formulation and performance analysis

PDF

Computational investigation of cholinergic modulation of the hippocampus and directional encoding of touch in somatosensory cortex

PDF

A double-layer multi-resolution classification model for decoding time-frequency features of hippocampal local field potential

PDF

Decoding memory from spatio-temporal patterns of neuronal spikes

PDF

Acquisition of human tissue elasticity properties using pressure sensors

PDF

Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations

PDF

Finite element model to understand the role of fingerprints in generating vibrations

PDF

Bidirectional neural interfaces for neuroprosthetics

PDF

Computational investigation of glutamatergic synaptic dynamics: role of ionotropic receptor distribution and astrocytic modulation of neuronal spike timing

Asset Metadata

Creator Girard, Christopher Benedict Charles (author)

Core Title Adaptive octree meshing of the extracellular space in neural simulations

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Biomedical Engineering

Degree Conferral Date 2023-08

Publication Date 07/26/2023

Defense Date 07/21/2023

Publisher University of Southern California. Libraries (digital)

Tag adaptive mesh,admittance method,compartmental model,finite-element,local field potential,OAI-PMH Harvest,volume conductor

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Song, Dong (committee chair), Bouteiller, Jean-Marie (committee member), Lazzi, Gianluca (committee member), Zhou, Qifa (committee member)

Creator Email cbcgirard@gmail.com,cgirard@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC113289506

Unique identifier UC113289506

Identifier etd-GirardChri-12152.pdf (filename)

Legacy Identifier etd-GirardChri-12152

Document Type Dissertation

Rights Girard, Christopher Benedict Charles

Internet Media Type application/pdf

Type texts

Source 20230726-usctheses-batch-1074 (batch), 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 author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA

Repository Email cisadmin@lib.usc.edu