Enabling Electric Field Model of Microscopically Realistic Brain

Enabling Electric Field Model of Microscopically Realistic Brain | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Enabling Electric Field Model of Microscopically Realistic Brain Zhen Qi , View ORCID Profile Gregory M. Noetscher , Alton Miles , Konstantin Weise , Thomas R. Knösche , Cameron R. Cadman , Alina R. Potashinsky , Kelu Liu , William A. Wartman , View ORCID Profile Guillermo Nunez Ponasso , Marom Bikson , Hanbing Lu , View ORCID Profile Zhi-De Deng , Aapo R. Nummenmaa , Sergey N. Makaroff doi: https://doi.org/10.1101/2024.04.04.588004 Zhen Qi 1 Department of Electrical and Computer Eng., Worcester Polytechnic Inst. , Worcester MA USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Gregory M. Noetscher 1 Department of Electrical and Computer Eng., Worcester Polytechnic Inst. , Worcester MA USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Gregory M. Noetscher For correspondence: gregn{at}wpi.edu Alton Miles 1 Department of Electrical and Computer Eng., Worcester Polytechnic Inst. , Worcester MA USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Konstantin Weise 2 Max Planck Inst. for Human Cognitive and Brain Sciences , Leipzig, Germany 3 Leipzig University of Applied Sciences (HTWK), Faculty of Engineering , Leipzig, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site Thomas R. Knösche 2 Max Planck Inst. for Human Cognitive and Brain Sciences , Leipzig, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site Cameron R. Cadman 1 Department of Electrical and Computer Eng., Worcester Polytechnic Inst. , Worcester MA USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Alina R. Potashinsky 1 Department of Electrical and Computer Eng., Worcester Polytechnic Inst. , Worcester MA USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Kelu Liu 1 Department of Electrical and Computer Eng., Worcester Polytechnic Inst. , Worcester MA USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site William A. Wartman 1 Department of Electrical and Computer Eng., Worcester Polytechnic Inst. , Worcester MA USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Guillermo Nunez Ponasso 1 Department of Electrical and Computer Eng., Worcester Polytechnic Inst. , Worcester MA USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Guillermo Nunez Ponasso Marom Bikson 4 Department of Biomedical Engineering, The City College of New York , New York NY USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Hanbing Lu 5 Neuroimaging Research Branch, National Institute on Drug Abuse, Intramural Research Program, National Institutes of Health , Baltimore MD USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Zhi-De Deng 6 Computational Neurostimulation Research Program, Noninvasive Neuromodulation Unit, Experimental Therapeutics and Pathophysiology Branch, National Institute of Mental Health, National Institutes of Health , Bethesda MD USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Zhi-De Deng Aapo R. Nummenmaa 7 Athinoula A. Martinos Ctr. for Biomedical Imaging, Massachusetts General Hospital , Charlestown MA USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Sergey N. Makaroff 1 Department of Electrical and Computer Eng., Worcester Polytechnic Inst. , Worcester MA USA 8 Department of Mathematical Sciences , Worcester Polytechnic Inst., Worcester MA USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Preview PDF Abstract Modeling brain stimulation at the microscopic scale may reveal new paradigms for various stimulation modalities. We present the largest map to date of extracellular electric field distributions within a layer L2/L3 mouse primary visual cortex brain sample. This was enabled by the automated analysis of serial section electron microscopy images with improved handling of image defects, covering a volume of 250 × 140 × 90 μm 3 . The map was obtained by applying a uniform brain stimulation electric field at three different polarizations and accurately computing microscopic field perturbations using the boundary element fast multipole method. We used the map to identify the effect of microscopic field perturbations on the activation thresholds of individual neurons. Previous relevant studies modeled a macroscopically homogeneous cortical volume. Our result shows that the microscopic field perturbations – an ‘electric field spatial noise’ with a mean value of zero – only modestly influence the macroscopically predicted stimulation field strengths necessary for neuronal activation. The thresholds do not change by more than 10% on average. Under the stated limitations and assumptions of our method, this result justifies the conventional theory of “invisible” neurons embedded in a macroscopic brain model for transcranial magnetic and transcranial electrical stimulation. However, our result is solely sample-specific and largely neglects the effect of the microcapillary network. Furthermore, we only considered the uniform impressed field and a single- pulse stimulation time course. Significance statement This study is arguably the first attempt to model brain stimulation at the microscopic scale, enabled by automated analysis of modern scanning electron microscopy images of the brain. It concentrates on modeling microscopic perturbations of the extracellular electric field caused by the physical cell structure and is applicable to any type of brain stimulation. Data availability statement Post-processed cell CAD models (383, stl format), microcapillary CAD models (34, stl format), post-processed neuron morphologies (267, swc format), extracellular electric field and potential distributions at different polarizations (267x3, MATLAB format), *.ses projects files for biophysical modeling with Neuron software (267x2, Neuron format), and computed neuron activating thresholds at different conditions (267x8, Excel tables, without the sample polarization correction from Section 2.8 ) are made available online through BossDB , a volumetric open-source database for 3D and 4D neuroscience data. 1. Introduction Advances in brain stimulation leverage models of electric fields within the brain for understanding the brain’s response to external stimuli provided by transcranial magnetic stimulation (TMS) [ 1 ], transcranial electrical stimulation (TES) [ 2 ], deep brain stimulation (DBS) [ 3 ], electroconvulsive therapy (ECT) [ 4 ], or similar excitation. However, accurately modeling these fields within the brain presents a great challenge due to the brain’s complex microstructure, including the intricate networks of neurons and blood vessels [ 5 ],[ 6 ],[ 7 ],[ 8 ],[ 9 ] all of which must be included into the model. The conventional approach of field modeling in brain stimulation [ 10 ],[ 11 ],[ 12 ],[ 13 ],[ 14 ],[ 15 ],[ 16 ] ignores all microscopic structures, with each region represented by macroscopically averaged properties (e.g., isotropic or anisotropic macroscopic conductivity). Only after the macroscopic electric field is evaluated with finite element method (FEM) [ 17 ], are isolated neurons “reinserted” to model their biophysical response. Ignoring microstructure is a recognized limitation but treated as a fait accompli because of known computational intractability of FEM for very complex geometries, except for models of isolated and simplified axonal and vascular geometries (cylinders, spheres, regular polygons, etc.), where the importance is demonstrated both theoretically and experimentally [ 18 ],[ 19 ],[ 20 ],[ 21 ],[ 22 ],[ 23 ],[ 24 ],[ 25 ],[ 26 ]. Over the last few years, several brain samples with realistic anatomical microstructure enabled by automated analysis of serial section electron microscopy images were made available. Those are an L2/L3 mouse primary visual cortex brain sample (250 ×140×90 μm 3 volume) or IARPA Phase I [ 5 ],[ 6 ],[ 7 ] with both membrane meshes and morphological reconstructions, a 1 mm^3 mouse primary visual cortex brain sample with membrane meshes and cortical column reconstruction or IARPA Phase III [ 7 ],[ 8 ], and a 1 mm^3 human brain sample (temporal cortex) with membrane meshes and some morphological reconstructions [ 9 ]. All samples are augmented with capillary blood networks (about 230 millimeters of blood vessels in [ 9 ]). Is it possible to electromagnetically model these samples to extract information useful for brain stimulation? The present study introduces arguably a first attempt to do so. In contrast to the common FEM, we utilize a matrix-free iterative boundary element fast multipole method (BEM- FMM) [ 27 ],[ 28 ],[ 29 ] with a new nested iterative solution based on physical separability of different cells. The BEM-FMM charge-based formulation [ 30 ],[ 31 ],[ 32 ],[ 33 ] could presently handle up to ∼200 million surface degrees of freedom simultaneously. It is readily applicable to computing electric charges induced on nonconducting cell membranes (under initial polarization) of realistic closely spaced neuronal cells and nonconducting (to the main order of approximation) walls of microcapillaries. These charges generate their own secondary electric field E s which alters the primary macroscopic brain stimulation field E i (due to a TMS coil for example) and thus changes expected activating threshold of the neurons. In this approximation, the extracellular solution becomes fully decoupled from the intracellular solution. For a finite brain sample subject to an impressed field E i , the extracellular solution is a unique solution of the exterior Neemann problem [ 37 ],[ 38 ]. Field variations on the cellular level result in an ‘electric field spatial noise’ with a mean value of zero. In contrast to temporal noise studies (cf. [ 39 ],[ 40 ]), such a noise is large and rather multiplicative since microscopic field perturbations are directly proportional to the local values of the impressed field E i itself. As a result, a map of distributions of the extracellular electric field within the entire brain sample has been created. The electric field computations alone required approximately 3 months of continuous CPU time, utilizing 7 modern terminal servers simultaneously as well as multiple nodes of a high-performance computing cluster. We used the map to identify microscopic perturbations of the extracellular electric field and their effect on the activating thresholds of individual neurons. Previous relevant studies modeled a macroscopically homogeneous cortical volume. To do so, we averaged the computed total extracellular field over membrane cross- sections and obtained a microscopically altered 1D activating function [ 41 ],[ 42 ],[ 43 ] at every cross-section of a neuronal process, in the approximation of the 1D cable equation. Thus, we did not model the volumetric intracellular space and did not employ the true bidomain approach [ 21 ],[ 44 ]. This is because the length-to-diameter ratio for an axonal arbor (most sensitive to activation) of a realistic cortical neuron maybe as high as 60,000 (based on data from [ 7 ],[ 8 ][ 9 ]), which facilitates the use of the 1D cable model. The extracellular electric field noise defined above alters the activating function, potentially causing local membrane depolarization and inward transmembrane current injection at its positive peaks, while nothing flows at the hyperpolarization peaks – the noise is being rectified. Therefore, our initial hypothesis was that microscopic extracellular noise substantially eases neuron activation as compared to microscopically homogeneous brain modeling. The question we would like to answer here is: by how much? To do so, standard biophysical analysis with Neuron software [ 45 ] was further carried out based on neuron morphologies precisely linked to the physical neurons and known ion channel models for rodents [ 46 ],[ 10 ],[ 11 ],[ 12 ]. The brain sample under study was considered as it stands i.e., no cloning, elongation, or scaling neuronal processes/capillaries of any kind are made. The study is organized as follows. Section 2 Materials and Methods describes assembly of the computational brain sample model (further detailed in Supplements A, B), describes computational workflow for the extracellular field (further verified in Supplement E) as well as the biophysical modeling setup (further detailed in Supplement C). Section 3 Results presents major computational findings for corrected activating thresholds with other results reported in Supplement C. Section 4 Discussion summarizes and explains practical findings of the study and states limitations of our approach. Section 5 Conclusions concludes the paper. 2. Materials and Methods 2.1 Brain sample – MICrONS L2/L3 dataset of mouse visual cortex The brain sample is a semi-automated reconstruction of the L2/3 P36 male mouse primary visual cortex volume of approximately 250×140×90 μm in size ( Fig. 1a-d ) from electron microscopic images with a resolution of 3.6×3.6×40 nm [ 5 ],[ 6 ],[ 7 ]. This reconstruction includes physical (not skeletonized) pyramidal and non-pyramidal neurons, astrocytes, microglia, oligodendrocytes and precursors, pericytes, vasculature, nuclei, mitochondria, and synapses [ 5 ],[ 7 ]. In total, 364 excitatory ( Fig. 1e ) and 32 inhibitory ( Fig. 1f ) neurons with cell bodies within the volume ( Fig. 1c ) were segmented, along with some microvasculature ( Fig. 1d ) [ 7 ]. For 301 neurons, the labeled neuronal centerlines (skeletons) in the common swc format were also made available [ 7 ]. This provides a unique opportunity to estimate the effect of the local field perturbations using biophysical analysis [ 45 ]. Download figure Open in new tab Fig. 1. Structure of L2L3 sample under study. a) Snapshot of relative position of a neuron (soma) vs. a nearby blood capillary (white area) within the sample [ 7 ], [ 5 ]. b) Snapshot of different neuronal processes labeled with different colors at a higher resolution. c) Segmented neuronal bodies of the entire sample. d) Centerlines/skeletons of individual neurons and segmented blood vessels in red (with an astrocyte at top right). The blood network of this sample was not the major target and is severely underdeveloped. e,f) Excitatory (pyramidal) and inhibitory (basket) neurons within the sample volume to scale. The axes differ from the mouse-centric coordinate system [ 5 ]. 2.2 Construction of computational cell meshes. Volumetric fraction of neurons Only the neurons themselves and the available vasculature (with a few astrocytes) – the largest perturbers of the impressed or primary electric field E i – were included into the present numerical analysis. Numerous neuroglia cells seen in Fig. 1a,b were not included in the analysis due to their small size, which made them impossible to model, even with advanced numerical simulation tools. The volume occupied by these microstructures (extracellular space) was considered approximately homogeneous. The volumetric fraction of the neuronal cells within the brain volume was estimated in Table S1 of Suppl. A as 30±5% by computing the volume of each neuron and taking into account neuronal mesh intersections (below). This is significantly less than the previously reported values on the order of 70% [ 48 ],[ 49 ],[ 50 ] (Suppl. A). One reason for it is that all neurons with somas outside of this relatively small sample were omitted and their dendritic/axonal processes, which are still within the sample, were ignored during the segmentation. Supplement A describes the post-processing procedure for the original sample and its results in detail. Based on this data, Table S1 of Supplement A lists major topological and biophysical parameter of the created computational testbed, including estimated sample anisotropy. Both the facet size of the triangular membrane discretization ( Fig. 2a ) for the electromagnetic analysis as well as the skeleton’s spatial resolution ( Fig. 2b,c ) for the biophysical analysis were kept at approximately 100 nm to resolve large yet very localized variations of the impressed field due to multiple nearfield interactions of the neuronal processes as well as multiple intersections of neuronal arbors of different neurons (below). Download figure Open in new tab Fig. 2. a) Final high-quality membrane surface triangular mesh for one neuron (#238) which is located approximately in the middle of the sample. The membrane mesh has 2 million facets in total and a 100 nm resolution. b,c) Healed centerlines of a neuron. d) Most common intersections observed – dendrite to dendrite and axon to dendrite intersections. Two different colors correspond to the two different neurons. e) Less common intersections – dendrite to soma and axon to axon, respectively. f,g) Mesh intersection resolution via removal of inner triangles. h) Result of this operation for an axon-to-dendrite intersection. Removed triangular facets are marked dark green. 2.3 Neuronal mesh intersections (synaptic interconnections) and their resolution The neuronal and vasculature surface meshes so constructed were found to heavily intersect with each other at multiple locations (cf. Fig. S2 of Suppl. A). Some of these intersections reflect true synaptic interconnections between different neurons with a list of approximately 2,000 synaptic clefts quantified in [ 6 ] for the same sample. At the same time, many intersections are likely due to the imprecise segmentation and subsequent, not entirely identical, remeshing. Fig. 2d illustrates the most common intersections – interconnections – including the dendrite to dendrite and axon to dendrite intersections. Two different colors correspond to the two different neurons. Fig. 2e illustrates the less common intersections – dendrite to soma and axon to axon, respectively. Additionally, Supplement A demonstrates an intersection map between neighboring excitatory neurons #238 and #290, respectively ( Fig. S2 ). Fig. S2 also presents a mesh intersection roster of excitatory neuron #238 with all other neurons of the sample. It is observed that neuron #238 intersects with almost all neurons of the sample. Similar results were observed for other larger neurons within the sample. To resolve the mesh intersections numerically, we used the method illustrated in Fig. 2f,g ,h. For two intersecting processes from Fig. 2f , all facets belonging to either process and having at least one node within another process (its surface mesh) have been removed. This operation does not exactly join the two meshes but leaves a narrow gap (one the order of triangle size or ∼100 nm) in between as seen in Fig. 2h . When the gap size is small compared to the intersection domain itself, which is achieved by a present fine mesh resolution, the gap has a minor influence on the results as demonstrated by our refined mesh computations (Supplement A). Note that this method removes coincident or interpenetrating membranes and joins the corresponding intracellular spaces. 2.4 Computing extracellular electric field and extracellular electric potential Our electromagnetic modeling setup only computes the extracellular electric field in a finite brain specimen. For the computation of the extracellular electric field, the approach developed in two relevant yet smaller-scale studies [ 32 ],[ 33 ] has been used. It assumes that the membrane is non- conducting (insulating) at the end of initial polarization [ 34 ],[ 55 ],[ 56 ], during the application of a relatively short (0.1-1 ms or so) brain stimulation pulse. The same condition applies to the walls of the blood capillaries. The condition of the insulating membrane makes it possible to fully decouple the extracellular solution in the form of a classic uniquely solved exterior problem with the simple Neumann (zero normal current density) boundary conditions at the outer surface of the membrane. In response to an impressed or primary electric field E i with the (local) potential φ i , induced surface charges will accumulate on the membrane surfaces. The secondary field of these charges, E s , will alter the extracellular electric field. The total field in the extracellular space will become E e = E i + E s and the total extracellular potential will become φ e = φ i + φ s . The goal of the electromagnetic part of this modeling pipeline is to find the secondary field and potential. We have applied a uniform electric field E i with the amplitude of 100 V/m to the entire sample for each independent field polarization ( x, y, z ). We assumed that the sample, which has the form of the non-orthogonal parallelepiped ( Fig. 1c ), is embedded into a homogeneous medium with the same extracellular conductivity. In this formulation, the exact conductivity value does not matter, and no box around the irregular sample is needed. However, this formulation still needs adjustment to match the impressed uniform electric field outside the sample with minimal global distortion. In other words, the finite microscopic sample must be embedded into a conventional microscopically homogeneous brain model with minimal distortion. There are different ways to accomplish this. One way is to increase extracellular conductivity within the finite specimen. This will lead to lowering the average extracellular electric field within the specimen. The method for determining this adjustment is based on modeling an artificial rectangular-cuboid sample that mimics the original specimen. This method is described in detail further in this section. 2.5 Application of BEM-FMM. Nested iterative algorithm The single-layer charge formulation for the exterior Neumann problem allows for direct application of the charge-based boundary element fast multipole method (BEM-FMM) developed previously [ 30 ],[ 31 ],[ 32 ],[ 33 ]. The BEM-FMM is ideal for this type of problem since it does not require construction of an explicit system matrix and does not require a volumetric mesh for the extremely complicated interior/exterior domains. For the present case, however, the iterative BEM-FMM was only capable of computing the field of ∼100 interacting neurons (with ca 0.2 billion facets) simultaneously. For larger samples, it required a prohibitively large RAM more than 1 TB and was experiencing slower convergence of the underlying GMRES algorithm [ 60 ]. Therefore, a nested iterative method has been developed and is illustrated in Fig. 3a . At initial step 1, each full neuron is solved in isolation and is not affected by charges deposited on other neurons. The membrane charge deposition in response to the impressed field E s (or the self-interaction) is found via the standard iterative BEM-FMM solution for every neuron in homogeneous space, with other neurons still no existing. Local charge density variations at the soma seen in Fig. 3a top row are solely due to the realistic irregular soma shapes, which are quite different from simplified spherical/ellipsoidal approximations. After that, all faces within other neurons have been removed ( Section 2.3 ). Download figure Open in new tab Fig. 3. a) Nested iterative solution method. At step 1, every neuron is considered independent. Membrane charge deposition in response to impressed or primary field E i ( E e i in the extracellular space) is found via an iterative BEM-FMM solution. At step 2, the primary field E e i is modified by the field of charges deposited on all other neurons. Membrane charge deposition for every neuron in response to this modified field is found next. The process further repeats itself until a convergence criterion is met. b) Convergence curves for the relative 2-norm surface charge density error. c) Convergence curves for the relative 2-norm error for the collinear electric field projected onto centerlines of neuronal processes. At step 2, the primary field E i ( E e i in the extracellular domain) is modified by the field of charges deposited on all other neurons found in step 1. The membrane charge deposition for every neuron is recomputed with this modified primary field. At the following steps, the process repeats itself until a convergence criterion is met. In Fig. 3 , the convergence criterion used required the relative 2-norm induced surface charge density error to be below 3% on average and the relative 2-norm centerline electric field error to be below 1% on average. The corresponding convergence curves ( Fig. 3b,c ) indicate a sub-percent error for the collinear electric field at the neuronal centerlines after 6 iterations. An additional test of solution accuracy has been performed by comparing with the global BEM-FMM solution for 11 neighbor neurons treated as a single network. This test is described in Supplement E. There, the relative average charge density error for every neuron versus the global iterative or ‘exact’ solution for that single network after six iterations was approximately 1%. Emphasize that the neuron-to-neuron interactions are taken into account at every step of the nested iterative solution except trivial first iteration ( Fig. 3a ). In its turn, the nested iterative algorithm used in this study provides a complete solution; it rapidly converges to the exact result when the number of iterations becomes sufficient (as shown in Supplement E). The BEM-FMM solution is already highly parallelized, so we did not attempt additional parallelization of the nested iterative loop. An outer loop selects a neuron under study. An inner loop retrieves the charges on all other neurons found previously and computes their fields on the facets of the neuron under study (adds these fields up). Next, the neuron under study is solved, i.e., its new charge distribution is obtained, etc. This method is potentially capable of modeling an unlimited number of cells, although CPU time now increases as O ( N 2 ) where N is the total number of cells. Electromagnetic computational results of this study required approximately 3 months of continuous CPU time when simultaneously utilizing multiple terminal servers and nodes of a high-performance computing cluster. Fig. 4 illustrates the final induced charge map for the dorsal-ventral ( y ) polarization. It demonstrates the induced surface charge density for neuron #321 and its 200 nearest neighbors. Fig. 4a,b is the induced surface charge density for neuron #321 in sagittal and coronal planes. Fig. 4a is the induced charge density when only neuron #321 is present (iteration 1 or self- interaction). Fig. 4b is the induced charge density when all neurons and microcapillaries of the sample are present (iteration 6). Fig. 4c,d is the induced surface charge density for 200 nearest neighbors of neuron #321. Fig. 4c gives the data for the complete network of 200 nearest neighbors. Fig. 4d is a zoomed in area (x33) of the network in Fig. 4c with some inter-synaptic interconnections clearly visible. We were unable to plot electromagnetic data for more than 200 neurons simultaneously due to the limitations of the graphical software used. Download figure Open in new tab Fig. 4. a,b) Induced surface charge density for neuron #321. a) Induced charge density when only neuron #321 is present (iteration 1 or self-interaction). b) Induced charge density when all neurons and microcapillaries of the sample are present (iteration 6). c,d) Induced surface charge density for 200 nearest neighbors of neuron #321. c) Complete network of 200 nearest neighbors. d) Zoomed in area (x33) of the network in c). We were unable to simultaneously plot electromagnetic data for more than 200 neurons due to the limitations of the graphical software used. 2.6 Average extracellular potential computed at centerlines/skeletons of neuronal processes The activating function of the one-dimensional cable equation [ 41 ],[ 42 ],[ 43 ] needs extracellular potential or the collinear extracellular electric field at the centerlines of neuronal processes in the form of infinitesimally thin, transparent ‘wires’. For finite-thickness neuronal processes, these values are extracellular fields averaged over the cross-section. In the extracellular domain ( and only there ), the field solution obtained previously is equivalent to the solution for the entirely non- conducting interior. Therefore, the averaged values are most accurately obtained when the field and the potential are directly computed at the centerlines using the same extracellular solution which simultaneously holds for the entirely non-conducting interior [ 32 ],[ 33 ]. Fig. 5 illustrates the present averaging method ( Fig. 5b ) and the resulting effect of membrane charge deposition on the collinear electric field projected onto the centerlines of neuronal processes for neuron #321 ( Fig. 5a,c ). The sagittal plane view is shown corresponding to Fig. 4a,b left. Fig. 5a is the electric field projection in homogeneous space, when the neuronal membrane is not present per se and the same unform impressed electric field of 100 V/m is applied along the y axis (from dorsal to ventral). Fig. 5c is the averaged collinear electric field projection when all neuron membranes are present, and the primary field is distorted by the field of the induced charges. One observes that perturbations of the collinear electric field are quite significant. Moreover, they are constructively or “in phase” accumulated in space so that a positive perturbation is added to the originally positive field projection and vice versa. Download figure Open in new tab Fig. 5. Effect of membrane charge deposition on the collinear electric field projected onto the centerlines of neuronal processes for neuron #321. The sagittal plane view is shown corresponding to Fig. 4a,b . a) Electric field projection in homogeneous space, when the neuronal arbor is not present and a uniform impressed electric field of 100 V/m is applied along the y axis (from dorsal to ventral). b) Averaging method over a cross-section. c) Electric field projection when all neurons are present, and the primary field is distorted by the field of all induced charges. An appealing result of this approach is that both the transmembrane potential and the double layer density at the end of initial polarization are obtained automatically. Once the extracellular potential φ e is known, we subtract its value from the constant resting potential value inside the cell and obtain the transmembrane potential φ m . Its product with the electrical permittivity of the membrane is the density of the double layer of charges across the cell membrane. 2.7 Biophysical modeling We implemented rodent neuron models developed in Refs. [ 46 ],[ 10 ],[ 11 ],[ 12 ] with 6 active channels for the soma, 8 active channels for the axons, 1 active channel for the basal dendrites, and 3 active channels for the apical dendrites initially developed and used by Aberra et al. [ 10 ]. The NEURON software [ 45 ] with fine spatial and temporal resolutions has been employed to study neural activation in response to the extracellular potential/electric field computed previously at the centerlines of the neuronal processes. The neuron morphology for 267 available neuron models (including soma, axons, apical and basal dendrites) was loaded into NEURON by its Import3D tool. The implementation for the most important axonal channels was additionally verified against in-house MATLAB-based biophysical software. Every neuron was modeled independently, following the approach of Refs. [ 10 ],[ 11 ],[ 12 ]. Intercellular communications were not yet considered. We used four stimulation pulse forms ( Fig. S3 of Supplement C): (i) a monophasic positive/negative rectangular pulse starting at t = 0 and with a duration of 0.1 ms and; (ii) an experimentally measured positive/negative biphasic TMS pulse with an overall duration of 0.3 ms. The pulse amplitude gradually increased until the neuron activating criterion (the transmembrane voltage rises above zero anywhere in the arbor) was met. After that, the amplitude was decreased with a smaller step to pass the activating threshold in the opposite direction. This automated search procedure – the bisection method – requires approximately 15 steps and about 15-20 min per neuron on average. The action potentials were partially initiated at the axonal terminations following the previous observations [ 10 ],[ 11 ],[ 12 ]. For some neurons – primarily where the neuronal processes have been substantially cut out – it was also the axon hillock. The overall simulation time for all neurons with x , y , and z polarizations, for all pulse forms, and for both primary and total extracellular fields/potentials did not exceed two CPU months. Other details of the biophysical modeling setup are described in Supplement C. The L2/L3 dataset does not provide information about axon myelination. However, the rodent axons are partially myelinated [ 10 ],[ 11 ],[ 12 ]. To account for this, we also modeled the case of artificial myelination. All axonal branches were assigned a myelin sheath with a thickness of 0.2 times the local radius. The nodes of Ranvier were placed at bifurcations, terminations, and along the segmental length, 20 µm apart from each other. 2.8 Correction of average extracellular field due to finite sample size Replacing a block of homogeneous tissue with the microscopic realistic specimen in a conventional macroscopic brain model must not change the solution in the surrounding homogeneous tissue. This in particular means that the extracellular conductivity of the heterogeneous specimen must be higher than the surrounding homogeneous conductivity, to enable nearly distortion free current flow through the specimen with non-conducting membranes, which are partially blocking the extracellular current flow. To do so, a box around the specimen is needed. For the present specimen (the non-orthogonal parallelepiped in Fig. 1c ), it is difficult to put it in a tightly enclosing rectangular box without substantial (up to 25%) reduction of the already small cellular volume, closing hundreds of holes, etc. Therefore, an artificial rectangular cuboid sample mimicking the original specimen was constructed as shown in Fig. 6 . Fig. 6a shows an equivalent straight fiber density distribution (computed in Supplement A) of the sample to be converted to the form of a filled rectangular cuboid. The fiber spacing there schematically illustrates weak sample anisotropy. Fig. 6b is an artificial computational equivalent of the specimen from Fig. 6a constructed using ∼500 finite cylinders in the form of a 3D lattice, with the overall mesh size of ca 2 million facets. The cylinders are spaced equally, but they have proportionally different radii to account for the moderate sample anisotropy from Fig. 6a . The cylinders are embedded into an enclosing (yellow) box (250×140×90 µm) whose conductivity contrast (interior vs. exterior conductivity) can be arbitrarily varied. The volumetric fill factor of the cylindrical lattice within the box is 0.3, which coincides with the fill factor for the original sample (Supplement A). The average extracellular field within the box is evaluated for 0.7 million uniformly distributed observation points. Cylinder walls (i.e., ‘membranes’) and their tips, which are protruding the box as I Fig. 6b , are assumed to be non-conducting. The field within the cylinders with the non-conducting walls (in the ‘intracellular’ space) is not computed. Download figure Open in new tab Fig. 6. a) Equivalent straight fiber density distribution of the sample converted to the form of a rectangular cuboid. The fiber spacing illustrates weak sample anisotropy. b) Artificial computational equivalent of the specimen from a) constructed using ∼500 finite cylinders in the form of a 3D lattice. The cylinders are spaced equally but have proportionally different radii to account for the moderate sample anisotropy from a). The cylinders are embedded into an enclosing (yellow) box (250×140×90 µm) whose conductivity contrast (interior vs. exterior conductivity) is varied. The volumetric fill factor of the cylindrical lattice within the box is 0.3, which coincides with the value for the original sample. Cylinder walls (membranes) and their tips, which are protruding the box, are assumed to be non-conducting. c-h) Magnitude (in V/m) and direction of the total extracellular electric field within and around the artificial sample for three different polarizations ( x , y , z ) of the impressed field E i with the value of 100 V/m. c,e,g) – Field distortion around the heterogeneous artificial sample when its extracellular conductivity is that of the surrounding homogeneous medium. d,f,h) – Field distortion around the heterogeneous artificial sample when its extracellular conductivity is increased with an attempt to achieve a minimum surrounding distortion. Fig. 6c-h shows the magnitude (in V/m) and direction of the total extracellular electric field within and around the artificial sample for three different polarizations ( x , y , z ) of the impressed field E i with a value of 100 V/m. Fig. 6c,e ,g illustrates the field distortion around the heterogeneous artificial sample when its extracellular conductivity is that of the surrounding homogeneous medium. Fig. 6d,f ,h provides the field distortion around the heterogeneous artificial sample when its extracellular conductivity is increased with an attempt to achieve a minimum overall distortion in the surrounding space. The external field distortion was estimated in three orthogonal planes shown in Fig. 6c-h . The extracellular conductivity was varied with the goal of minimizing the distortion over the outer plane area for all three planes simultaneously. An optimal extracellular conductivity increase by a factor of 1.7 was found using brute force optimization. Note that this conductivity increase within the sample appears to be close to that predicted by the Maxwell equation for dilute cell suspensions [ 57 ],[ 58 ],[ 59 ]. The correction factor is the ratio of the average extracellular field magnitude within the sample with the original and increased conductivity, respectively. This ratio was found to be equal to 1.08 for the x -polarization, 1.14 for the y -polarization, and 1.23 for the z -polarization. The activating thresholds found in the original specimen were multiplied by those values to account for sample matching to the surrounding homogeneous space. 3. Results 3.1 Neuronal activating thresholds within the sample vs homogeneous medium We chose the first 90 neurons with the longest non-cut axonal arbors exceeding 0.5 mm (Supplement A). The neurons have been renumbered in descending order based on the length of their axonal arbors. This operation eliminates, to the best of our abilities, neurons with severely truncated axonal arbors (at the sample boundaries) from the subsequent analysis, which would otherwise yield non-physical results. Fig. 7 shows the computed activating thresholds within the sample matched to the surrounding homogeneous medium compared to those in the microscopically homogeneous medium. Here and in the following, we only consider the 90 “best” neurons. Download figure Open in new tab Fig. 7. a) Equivalent straight fiber density distribution of the sample characterizing its anisotropy. b) Activating threshold reduction factor for the first 90 neurons with the longest axonal arbors (sorted) for three polarizations of the primary electric field: medial-lateral (red), dorsal-ventral (green), and anterior-posterior (blue). c) Absolute activating thresholds for the same 90 neurons in homogeneous space. d) Absolute activating thresholds for the same 90 neurons within the sample. The units in c,d) are given in terms of the base field strength of 100 V/m; one unit corresponds to 100 V/m. These results are given for the case of the unmyelinated axonal arbor and a negative rectangular activation monophasic pulse with the duration of 0.1 ms. Results for all other cases are given in Supplement C. Fig. 7b gives a neuron-by-neuron threshold reduction factor. The reduction factor is defined here as the ratio of the activating threshold of a physical neuron within the sample (with the total extracellular field E e and the total potential φ e ) to the activating threshold of the neuron skeleton in the homogeneous medium (subject to the field E i and the potential φ i ). The threshold reduction results in Fig. 7b are given for three polarizations of the primary electric field: medial- lateral (red), dorsal-ventral (green), and anterior-posterior (blue) and for the 0.1 ms long rectangular negative stimulating pulse. Fig. 7c gives the corresponding absolute activating thresholds for 90 neuron skeletons in the homogeneous space. Fig. 7d gives the absolute activating thresholds for the same 90 physical neurons within the sample. The units in Fig. 7c,d are given in terms of the field strength of 100 V/m; one unit corresponds to 100 V/m. Figs. S4-S10 of Supplement C report the same data for the remaining pulse forms and for myelinated axons. The last figure there was computed for the positive biphasic 0.3 ms long pulse. The data in Supplement C are conceptually similar to those shown in Fig. 7 . Average threshold reduction factors for all pulse forms and polarizations of E i are summarized in Table 1 . View this table: View inline View popup Download powerpoint Table 1. Averaged relative threshold reduction factor of the neuronal activating thresholds for different polarizations of the impressed (primary) electric field and different pulse forms (Supplement C) between a neuronal skeleton in a homogeneous medium and the physical neuron in the brain sample, respectively. The first 90 neurons with the longest (uncut) axonal arbors have been processed. Since data for myelinated axons are quite similar (Supplement C), only one such dataset has been included. Correction for the average field was found by modeling an artificial rectangular cuboid of the same effective size and fill factor including anisotropy from Fig. 6b ( Section 2.8 ). 3.2 Deviation of activating thresholds for different polarizations of primary field Table 2 summarizes normalized standard deviations of the average neuronal activating thresholds for three different polarizations of the primary electric field. They are computed as std([ X Y Z ])/mean([ X Y Z ]) where X Y Z are the average thresholds for three respective polarizations. The first 90 neurons with the longest axonal arbor have again been processed. From Table 2 and Fig. 7 , one observes that the normalized standard deviations between average thresholds for different directions of the applied primary field are always smaller for the sample as compared to individual neuron skeletons in homogeneous space. Cell membrane interactions apparently tend to make brain stimulation more ‘equal’ in all directions. View this table: View inline View popup Download powerpoint Table 2. Normalized standard deviation of the neuronal activating thresholds for three different polarizations of the primary electric field (from Fig. 7 and Supplement C). The first 90 neurons with the longest (largely uncut) axonal arbors have been processed. 3.3 Strength-duration results Fig. 8a illustrates a weak dependency of the threshold reduction factor on the duration of the applied stimulation pulse. A case in point is neuron #348 within the sample when the negative rectangular pulse is applied. At the same time, the absolute values of the activating threshold indeed follow the classic strength duration dependency [ 61 ] with a high degree of accuracy. This is shown in the inset in Fig. 8a ; note that a log-log scale has been used there. Download figure Open in new tab Fig. 8. a) Activating threshold reduction factors as functions of rectangular negative pulse duration. The inset shows the corresponding strength-duration curves for neuron #348 within the brain sample. Inset units in are in terms of the initial field strength of 100 V/m; one unit corresponds to 100 V/m. b,c,d) Activating threshold reduction factors as functions of the number of neighboring cells for three representative excitatory neurons with different lengths of the axonal arbor (for the rectangular negative pulse) without the sample polarization correction. All three neurons belong to the set of 90 neurons with the longest axonal arbor. The x-axis is the corresponding volumetric fraction of neuronal mass computed separately for each neuron when sequentially adding its neighboring cells. Red lines show the expected fraction of 70% [ 50 ]; dashed lines show expected extrapolation results. Vertical black lines on the right and the corresponding larger circles are the results with the sample polarization correction for volumetric concentration of 0.5. 3.4 Dynamics of threshold reduction factor as a function of a number of surrounding neurons Fig. 8b,c ,d illustrates the threshold reduction factor dynamics as a function of the number of neighboring cells included in the computations for three representative neurons located close enough to the sample center. For each of three representative neuron cells (#238, 321, 348, all having different lengths of axonal arbor), we start adding its nearest neighbors (0, 5, 10 neighbors, etc.) from the entire brain sample and then repeat both electromagnetic and biophysical computations to find new activating thresholds. The x -axis in Fig. 8b,c ,d is the corresponding relative volumetric fraction of neuronal mass computed using known volumes of each cell. The first point of every curve is the threshold reduction factor without neighbors, which is due to self-interaction of the neuron itself; the last point is the result for complete brain sample. These curves do not yet include the sample polarization correction, which was very difficult to obtain for the scattered neurons. 3.5 Extrapolation of threshold reduction factor to more realistic neuronal density The volumetric fraction of the neuronal cells within the volume was estimated in Table S1 (Supplement A) as ∼30%. This is less than the known value of ∼60-70% [ 48 ],[ 49 ],[ 50 ] in rats and mice since all neurons with somas outside the brain sample were omitted and their dendritic/axonal processes, which are still within the sample, were ignored. The threshold reduction factor curves in Fig. 8b,c ,d could in principle be extrapolated (as shown by dashed lines in Fig. 8b,c ,d) to estimate the threshold reduction factor for the likely more realistic volumetric neuronal fill factor of 70% (shown by a red vertical line in Fig. 8b,c ,d). To estimate the effect of the sample polarization correction, we again computed the expected extracellular field decrease for one volumetric concentration of 0.5 following the method of Section 2.8 . In Fig. 8b,c ,d, the vertical black lines on the right and the larger colored circles are the final results with the sample polarization correction. 4. Discussion 4.1. Methodology The present study is a method to electromagnetically model a realistic microscopic brain sample in application to brain stimulation. Our target region is a relatively small L2/L3 mouse primary visual cortex brain sample (250 ×140×90 μm 3 volume) or IARPA Phase I [ 5 ],[ 6 ],[ 7 ] with less than 400 tightly packed neurons, with much of the neuronal arbor missing due to excluding cells with somas located outside the volume under study, and with an underdeveloped microvasculature. We show that the boundary element fast multipole method or BEM-FMM [ 27 ],[ 28 ],[ 29 ] – when using its charge-based formulation [ 30 ],[ 31 ],[ 32 ],[ 33 ] – is an appropriate tool for the combined electromagnetic and biophysical analysis of not only this sample, but also potentially more accurate modern detailed brain samples [ 7 ],[ 8 ],[ 9 ] with the resolution of at least 100 nm under certain assumptions. These assumptions include: Only major scatterers or perturbers of the applied electric field – neurons themselves and microcapillaries – are included into consideration. The remaining cells are combined into one extracellular space. This assumption may be revised with regard to astrocytes and oligodendrocytes [ 9 ]. Cell membranes are assumed to be nonconducting at the end of the initial polarization period for a duration of several microseconds [ 34 ]. This assumption remains valid until the end of a stimulation pulse (with a duration of approximately 0.1-1 ms). In this time interval, we compute the activating function via the one dimensional cable equation [ 41 ],[ 42 ],[ 43 ] by averaging the total collinear extracellular field over a cross-section of a neuronal process. The neuronal cells then experience a change in their physiological state with the characteristic time scale of several milliseconds or longer [ 34 ]. At this point, we are again following the one-dimensional cable equation of the standard biophysical theory but with the activating function equal to zero. Our modeling method is not limited to the present case of a uniform applied brain stimulation field, which is more common for transcranial magnetic and electrical stimulation. It can equally well target any type of brain stimulation, including biophysical aspects of deep brain stimulation [ 66 ],[ 67 ],[ 68 ] and intracortical microsimulation/brain-computer interfaces (cf., for example [ 69 ],[ 70 ],[ 71 ]). 4.2. Lowering activating thresholds? Our initial hypothesis was that microscopic extracellular electric field noise might substantially ease neuron activation compared to macroscopically homogeneous brain modeling. However, this hypothesis was generally not confirmed by the present modeling results. More specifically: - The average activating thresholds remain nearly the same for electric field polarizations aligned along the two longest sample directions – medial-lateral and dorsal-ventral (see Table 1 ). - Some minor variations (reduction factors) observed there are primarily attributed to lowering of activating thresholds for neurons with already high absolute activating thresholds. Such neurons do not significantly contribute to recruitment. On the other hand, the neurons with the lowest activating thresholds – the most sensitive neurons (the ones of real interest) – are mostly unaffected (cf. Fig. 7 and the corresponding figures of Supplement C). - Although the average activating threshold for the anterior-posterior direction of the applied field does become lower, it is again significantly affected by lowering activating thresholds for neurons with already high absolute activating thresholds ( Fig. 7 and the corresponding figures of Supplement C). Furthermore, this is the shortest sample dimension. There, the absolute activating thresholds are high and the effect of sample polarization might not be entirely undone. - Extrapolation to a more realistic neuronal mass yields a similar result ( Fig. 8 ). Still, the fiber in the anterior-posterior direction is subject to the largest number of intersections with perpendicular fibers per unit length, as evidenced by Fig. 7a , and is proportional to the fiber density in the two other directions, i.e. to 1.0 e 5 × 1.3 e 5. For fibers in the dorsal-ventral direction, the number of such intersections per unit length in Fig. 7a is medium; it is again proportional to the fiber density in two remaining directions, i.e. to 0.9 e 5 × 1.3 e 5. Finally, for fibers in the medial- lateral direction, the number of such intersections per unit length in Fig. 7a is the lowest and is proportional to the fiber density in two other directions, i.e. to 0.9 e 5 × 1.0 e 5. These three numbers of intersections per unit length relate as 1.0, 0.90, and 0.70 and correlate well (Pearson correlation coefficient of 0.999) with the activating threshold reduction factors from Fig. 7b . Therefore, further work with larger regularly shaped samples may be necessary to prove or disprove this effect and estimate its strength. 4.3. When might field perturbations be important? Consider, for example, a simple numerical microsimulation experiment with a Hodgkin-Huxley axon originally studied by Rattay [ 43 ]. A cathodic stimulation of the 10 mm long axon by a point electrode located at the origin is shown in Fig. 9a . A monophasic electric pulse ( T = 0.1 ms) is applied. Field perturbers are 2D cylinders with non-conducting walls located perpendicular to the figure plane and with a large radius of 200 µm. Charge deposition on the cylinder walls perturbs the electrode field from Fig. 9a as shown in Fig. 9b where the total field is displayed. As a result, an electric field ‘spatial noise’ with a mean value of zero appears for the extracellular collinear electric field along the axon, as shown in Fig. 9c in blue. The noise amplitude is large; it can approach or even exceed the applied field itself when the scatterers are in close proximity to the axon (cf. basic solutions for DC current conduction, for example from [ 47 ]). Download figure Open in new tab Fig. 9. Concept of the electric field spatial noise. a) Straight HH axon from Ref. [ 43 ] and the electric field of cathodic stimulation. b) The same case as in a) but with artificial field perturbers – infinite cylinders. c,d) Collinear electric field and extracellular potential with (blue) and without (red) field perturbations. In both cases, the electrode current is 0.94 mA. e-n) Resulting evolution of transmembrane potential and transmembrane ionic current density for both cases at different times. Red curves are those without field perturbations; blue curves correspond to the perturbed extracellular field. When t < T ( Fig. 9e,f ), the transmembrane potential closely follows the activating function (the derivative of the collinear electric field) with sharp positive peaks of the latter indicating membrane depolarization [ 43 ]. The inward (negative in Fig. 9 ) ionic current flows within the depolarization peaks ( Fig. 1f ) while nothing flows within the hyperpolarization peaks – the noise is being ‘rectified’. This process is strongly nonlinear – cf. Fig. 9f,h . There, one local inward current injection starts to dominate other inward injections while all outward ionic currents are still negligibly small. Despite the following low pass filtering as in Fig. 9i,j , this local current injection appears to be strong enough to support further action potential generation ( Fig. 9k,l ,m,n) when the stimulation cathodic current is only 0.91 mA. This contrasts with the effect of the unperturbed field from Fig. 9a which would generate an action potential if the cathode current were 1.31 mA or higher. Thus, the threshold reduction factor due to the spatial noise is 0.91/1.31=0.69 in the present case. When the cylinder radius in Fig. 9b is halved (becomes 100 µm), the threshold reduction factor increases to 0.84. When the radius is further halved and is reduced to 50 µm, the threshold factor becomes 0.98 and then quickly approaches one when this process continues. Similar results hold for the mouse axon. It appears that strong yet very spatially narrow field perturbations caused by smaller cylinders do not have a significant effect, as they diffuse quickly along the cable with no lasting impact. This observation is directly evident during the simulations and is valid at the terminals as well. According to this example, to have an effect, the inward current injections must be not only strong but also sufficiently extended in space. The product of strength times the spatial length appears to be crucial. Proper lengths are observed in the case of blood vessels or other “unexpected” scatterers, such as possibly a combination of 2-3 nearby somas, a large astrocyte or oligodendrocyte, a favorable combination of several smaller scatterers, etc. According to the present simple example, the reduction effect becomes apparent only for relatively large scatterer sizes, on the order of ∼0.1λ where λ is the axonal space constant. This constant is equal to 1.5 mm for the Hodgkin-Huxley axon and to 0.6 mm for the mouse axon studied here. A similar observation holds for a finite axon (either Hodgkin-Huxley [ 43 ] or mouse [ 10 ],[ 11 ],[ 12 ]) in a uniform field subject to a pair of cylindrical scatterers at its termination or to a bandlimited electric field noise applied to an entire axon. The simple MATLAB code used to create this example is described in Supplement F and made available online [ 72 ]. 4.4 Limitations of this study The small size of the L2/L3 brain sample leads to cutting out a large portion of the neuronal processes, most notably axon collaterals. As an example, the average length of the existing axonal arbor per neuron is ∼0.5 mm for the present L2/L3 sample, ∼2.5 mm (with the maximum length of 36 mm!) for the 1mm^3 sample from Allen Inst. [ 7 ],[ 8 ], ∼6 mm for Markram’s et al elongated data in L2/L3, and ∼7 mm for Markram’s et al data in L5 [ 46 ]. Our computations show that the absolute activating threshold values are severely affected by the total axonal length. Additionally, the arguably most sensitive excitatory neurons of L5 [ 11 ] are not included into the sample. This is perhaps why the absolute threshold values in Fig. 7 and Figs. S4-S10 (Supplement C) overestimate the experimentally measured value of ∼90 V/m in mouse primary motor cortex (Supplement D) very significantly. Even in the most favorable case of myelinated axons, there are only nine neurons with the activating threshold below 200 V/m ( Fig. S8 , Supplement C) and none of them has the activating threshold below 100 V/m. Furthermore, the blood network of the L2/L3 sample was not the major target of IARPA Phase I study. It is therefore severely underdeveloped (cf. Fig. 1d ) compared to more recent and accurate results [ 7 ],[ 8 ],[ 9 ]. Our biophysical modeling treats all neuronal cells independently, following previous work [ 10 ],[ 11 ],[ 12 ] to demonstrate the impact of realistic microscopic structures on direct neuron activation threshold. On the other hand, the L2/L3 brain sample is already supplied with a list of ∼2,000 synapses between 334 excitatory pyramidal cells [ 7 ] obtained in [ 6 ] allowing future extension of our approach to simulating network/synaptic modulation [ 64 ]. Furthermore, we only considered the uniform impressed stimulation field and a single-pulse stimulation time course. Though used ubiquitously, there are inherent limitations of the one-dimensional cable equation, especially close to the large somas. Modern bidomain models (cf. [ 44 ]) could possibly resolve this issue. The next logical yet extremely challenging step is to analyze the much larger 1mm^3 MICrONS mouse brain sample, which includes a considerably better-developed axonal arbor and ∼75,000 neuronal cells, and another very recently published human brain sample [ 9 ]. This would require further improvements and modifications of the computational method. 5. Conclusions A method is proposed to accurately model perturbations of an impressed extracellular electric field within a microscopically realistic brain, with many tightly spaced neuronal cells. This impressed field could be due either to brain stimulation or to a firing cluster of neurons themselves. For the first time, we have been able to create a map of extracellular electric field distributions within the relatively sizable fully anatomical L2/L3 mouse brain sample (250 × 140 × 90 μm³) containing ∼400 detailed neuronal cells and some blood microcapillaries with the resolution of 100 nm. Electric field modeling was further coupled with 1D biophysical modeling, which is based on extracellular field averaging, neuron morphologies precisely linked to the physical neurons and known ion channel models for rodents. This combination was applied to predict changes in neuron activating thresholds caused by microenvironment – the electric field spatial noise with a mean value of zero on the cellular level. Under stated assumptions and limitations, our results show that the microscopic field perturbations – the electric field spatial noise caused by the neuronal structure – only modestly influence the predicted stimulation field strengths necessary for neuronal activation. This result supports the conventional theory of “invisible” neurons embedded in a macroscopic brain model [ 10 ],[ 11 ],[ 12 ], for transcranial magnetic and transcranial electrical stimulation. Although large in magnitude, the microscopic activating function perturbations are predominantly limited to very narrow spatial regions (e.g., dendritic diameters) so that the local inward current injections diffuse rather quickly instead of triggering substantial changes in the physiological state. However, our result is sample-specific and subject to a number of limitations and restrictions as stated above. Supplement A Post-processing segmentation and skeletonization of L2/L3 sample A1. Preparing neuronal surface meshes for electromagnetic modeling All 396 neurons were initially numbered following the y-coordinate (dorsal-to-ventral direction) of their center of mass. The lowest numbers correspond to the deepest neurons where the bulk of the axonal arbor was cut out (see below). Neurons with the highest numbers have the (relatively) well developed axonal arbor ( Fig. 1 of the main text). Initial computational surface triangular meshes for the neurons downloaded from the MICrONS database [ 1 ] were not entirely suitable for accurate numerical analysis. Automated CGAL alpha wrapping C++ library [ 2 ],[ 3 ] was used for all neurons to re-mesh the data and convert them to the strictly 2 manifold high-quality (min triangle quality is 0.1) high-resolution individual surface objects. A relative α (vs 1e-3*longest diagonal; parameter α limits the maximum triangle size) of 1000 and a relative offset (vs 1e-3*longest diagonal, parameter offset limits the tightness of the result around the input mesh) of 4000 were used in the CGAL program. We were unable to properly heal 13 cells. These cells were eliminated from the analysis which reduced the total cell number from 396 to 383. Fig. 2 of the main text illustrates the final membrane mesh for neuron #238. A typical membrane mesh has approximately 1-2 million triangles on average, with an average facet size kept at 100 nm. This is necessary to keep a proper spatial accuracy for intersection resolution ( Section 2.3 of the main text) and also resolve microscopic, localized variations of the impressed electric field. The total number of facets of the sample approaches 0.5 billion. The average facet size is 100 nm. An identical procedure also was applied to the microcapillaries (cf. Fig. 1d of the main text). A2. Volumetric density of neuronal cells computed from surface meshes and its comparison with literature values The volumetric fraction of all neuronal cells within the SEM (scanning electron microscopy) volume was estimated as 30%±5% by computing the volume of each neuron using the CGAL tools and also taking into account neuronal mesh intersections (Section 2.3 of the main text). This is less than the known estimates in rats and mice [ 4 ],[ 5 ],[ 6 ] because all neurons with somas outside the brain sample were omitted and their dendritic/axonal processes, which are still within the sample, were ignored. In general, the brain tissue is composed of extracellular space and intracellular space, i.e. neurons and glia cells (intracellular space) do not take up all the volume of brain tissue. Oxygen molecules and nutrients (e.g. glucose molecules) diffuse from capillary through extracellular space to neurons and glia, and are taken up through transporters. According to [ 4 ], neuronal elements occupy about 74% and glial elements about 8% of the total volume in rats. According to [ 5 ] (rats), the extracellular space takes up about 20% of total cortex volume in adult rats. According to [ 6 ] (Table 5), the fraction of neuronal mass in mice cerebral cortex is 70%. Based on these data and also on some data from Ref. [ 7 ], we estimate the realistic volume fraction of neuronal mass in mice neocortex as approximately 70%. A3. Healing and refining centerlines of neuronal processes – neuron skeletons The original neuronal skeletons or morphological reconstructions (301 in total) from the MICrONS database [ 1 ] didn’t remain entirely in the center of the corresponding neuronal processes represented by the relevant surface mesh. They also had a rather sparse discretization. The skeleton branches or centerlines have been therefore refined, smoothed, and corrected via an automated in-house routine that finds every point’s distance from the enclosing membrane walls in the up, down, left, right, in, and out directions and averages them for every cross-section of the relevant triangular mesh. Fig. 2b,c of the main text shows the centerlines of the neuronal processes healed in such a way. We attempted to maintain a nearly uniform spatial resolution of 100 nm for the centerline discretization(observation) points, which should be sufficient for catching up fast field variations over the length on the order of 1 µm or so. We were unable to properly heal 34 centerlines. These centerlines were eliminated from the subsequent analysis which reduced their total number from 301 to 267. We also estimated the average axonal length per neuron, with the maximum length of 8.3 mm for inhibitory neuron #191 and the minimum length of 0 mm (for 39 other neurons in total). The latter case is non-physical; it is a deficiency of the centerline extraction method. A4. Estimating sample anisotropy from neuron morphology To estimate the degree of anisotropy, we introduced an equivalent straight–fiber density measure. It is obtained when centerlines for all neuronal processes including both axons and dendrites are projected onto three Cartesian axes ( Fig. S1b ,c,d below). This density was estimated by finding the total projection length and then dividing this value by the product of three sample dimensions (the sample volume). The following values of the equivalent straight fiber density have been obtained: 1.3 e 5 mm −2 for the medial-lateral or x direction, 1.0 e 5 mm −2 for the dorsal-ventral or y direction and 0.9 e 5 mm −2 for the anterior-posterior or z direction. A5. Selecting proper neurons with the longest retained axonal arbors for biophysical modeling The lengths of all segments of the neuronal processes (axons and dendrites) were computed using the available *.swc files from the MICrONS database. Fig. S1a shows the length of axonal arbor per neuron, now sorted in descending order. It is seen that neurons with numbers greater than 90-100 have a vanishingly short axonal arbor, which was cut out of this sample almost entirely. Therefore, they were eliminated from the final results of biophysical analysis (otherwise non-physical, large values of the activating threshold could be obtained) but were still kept for the electromagnetic analysis. Thus, only the neurons with the overall axonal length of 0.5 mm or greater (with the sorted numbers from 1 to 90, Fig. S1a below) were included into the relevant biophysical modeling outcomes – see Fig. 6 of the main text, Figs. S4 through S7 of Supplement C, and also Tables 1 and 2 of the main text. Fig. S1b additionally shows the total length of axonal arbor per neuron projected onto three Cartesian axes (x, y, z). The neurons have been sorted as in Fig. S1a ). Similarly, Fig. S1c demonstrates projected (onto three Cartesian axes) length of dendritic arbor per neuron. The neurons have again been sorted as in Fig. S1a . Finally, Fig. S1d demonstrates the projected length of all neuronal processes (axonal plus dendritic) per neuron. The neurons have been sorted as in Fig. S1a . Again, one clearly observes in Fig. S1b that the neurons with numbers greater than approximately 90 have vanishingly small axonal arbor, which was cut out of the brain sample almost entirely. Download figure Open in new tab Fig. S1. a) Sorted (in descending order) length of axonal arbor per neuron. Neurons with numbers greater than approximately 90 have vanishingly small axonal arbor, which was cut out of the sample almost entirely. b) Projected (onto three Cartesian axes) length of axonal arbor per neuron. The neurons have been sorted as in a). c) Projected (onto three Cartesian axes) length of dendritic arbor per neuron. The neurons have been sorted as in a). c,d) Projected (onto three Cartesian axes) length of all neuronal processes (axonal plus dendritic) for all neurons with the healed morphological reconstructions. The neurons have been sorted as in a). A6. Summary of post-processing parameters Table S1 of this Supplement lists major topological and biophysical parameters of the created computational testbed (electromagnetic and biophysical) including estimated sample anisotropy. Emphasize again that both the facet size for membrane discretization and the skeleton’s spatial resolution for the biophysical analysis were kept at approximately 100 nm to resolve large yet very localized variations of the impressed field due to multiple nearfield interactions of the neuronal processes as well as to resolve multiple intersections of neuronal arbors of different neurons ( Fig. 2 of the main text). View this table: View inline View popup Download powerpoint Table S1. Post-processing parameters of the L2/L3 mouse brain sample.383 Neurons were kept for electromagnetic analysis. Data for 90 neurons with the longest axonal arbor were included in the final biophysical outcomes. Supplement B Membrane mesh intersection maps B1. Mesh intersection maps Fig. S2a demonstrates a typical surface mesh intersection map. Two neighboring excitatory neurons – #238 and #290 – have been chosen. Download figure Open in new tab Fig. S2. a) Surface mesh intersection map between neighboring excitatory neurons #238 and #290, respectively. Intersection domains are marked by magenta color. b) Mesh intersection roster of neuron #238 with all other neurons of the sample. The y axis gives the relative number of facets (% versus the total number of faces for neuron #238) with at least one facet node contained within another neuronal shell. The x axis is the neuron number. In Fig. S2a , the mesh intersection domains are marked by magenta color. There are approximately 30 such independent intersection domains. Fig. S2b also presents a mesh intersection roster of excitatory neuron #238 with all other neurons of the sample. It is observed in Fig. S2b that neuron #238 intersects with almost all neurons of the sample. Similar, and even more dramatic results were observed for other larger neurons of the sample. To identify the mesh intersections for large yet manifold meshes, we used vectorized function inpolyhedron written by Sven Holcombe in MATLAB [ 8 ]. Still, this has been a very time-consuming operation. Supplement C Biophysical modeling C1. Summary of NEURON modeling parameters Standard NEURON software [ 9 ] has been used to study neural activation in response to the impressed extracellular potential/electric field computed previously at the centerlines of the neuronal processes. The post-processed neuron morphology (described in Supplement A) for 267 available neuron models (in swc format, including soma, axons, apical and basal dendrites) was loaded into NEURON by its Import3D tool. Following references [ 10 ],[ 11 ],[ 12 ] we applied the existing ion channel models for rodents (rats) which have been used starting with Aberra et al. [ 10 ]. We used 6 active channels for the soma, 8 active channels for the axons, 1 active channel for the basal dendrites, and 3 active channels for the apical dendrites. The implementation for the most important axonal channels was tested and verified against in-house MATLAB-based software. The centerlines of the neuronal processes were discretized in space with the step d λ = L /[ n seg λ w ( f = 100 Hz )], where L is the section length, n seg is the number of segments in this section, and is the length constant, with D being the neuron diameter; f is a frequency that is high enough for transmembrane current to be primarily capacitive, but is still within the frequency range relevant to neuronal function. Here, R a is the resistivity of axoplasm and c m is the membrane capacity per unit area. In the computations, we tested several d λ values between d λ = 0.01 and d λ = 0.0005, respectively, to address expected rapid field variations at microscale. We compared the results for a representative number (6 in total) of large neurons with many local extracellular field variations. The difference in the activating thresholds was found to be negligibly small between the two extreme values. Therefore, the value d λ = 0.01 or 0.005 has been routinely used everywhere after test runs. The time step size was chosen as 1 μs which satisfies the stability condition. C2. Pulse forms of the stimulation field As the simple base form, we used the monophasic rectangular pulse form for the primary or impressed field starting at t = 0 and with a duration of 0.1 ms – cf. Fig. S3a ,b). The positive and negative pulses (positive and negative electric field directions) have been modeled separately. Along with this, we used a realistic biphasic TMS (transcranial magnetic stimulation) pulse form ( Fig. S3c ,d) from a B65 coil of MagVenture, Denmark measured in the Noninvasive Neuromodulation Unit, Experimental Therapeutics & Pathophysiology Branch, NIMH. Both positive and negative pulse forms have been investigated. Some computations have also been done for a monophasic TMS pulse. The pulse magnitude was gradually increased until the neuron activating criterion (the transmembrane voltage goes above zero anywhere in the arbor [ 10 ]) was met. After that, the amplitude was decreased with a smaller step to pass the activating threshold in the opposite direction. This automated search procedure – the bisection method – requires approximately 15 steps and about 15-25 min per neuron on average. The total simulation time for all neurons with x, y, and z polarizations for both primary and total extracellular fields/potentials did not exceed two weeks per polarization. Download figure Open in new tab Fig. S3. Four pulse forms for the primary impressed field used for biophysical modeling. C3. Reduction of neuronal activating thresholds for different pulse forms and myelinated axons Figs. S4 through S10 given below replicate Fig. 7 of the main text for different pulse forms and with considering axonal myelination. They have been used to fill out Tables 1 and 2 of the main text. Download figure Open in new tab Fig. S4. a) Equivalent straight fiber density distribution of the sample characterizing its anisotropy. b) Activating threshold reduction factor for the first 90 neurons with the longest axonal arbor (sorted) for three polarizations of the primary electric field: medial-lateral (red), dorsal-ventral (green), and anterior-posterior (blue). c) Absolute activating thresholds for the same 90 neurons in homogeneous space. d) Absolute activating thresholds for the same 90 neurons within the sample. The units in c,d) are given in terms of the base field strength of 100 V/m; one unit corresponds to 100 V/m. These results are given for the case of the unmyelinated axonal arbor and a positive rectangular activation monophasic pulse with the duration of 0.1 ms. Download figure Open in new tab Fig. S5. a) Equivalent straight fiber density distribution of the sample characterizing its anisotropy. b) Activating threshold reduction factor for the first 90 neurons with the longest axonal arbor (sorted) for three polarizations of the primary electric field: medial-lateral (red), dorsal-ventral (green), and anterior-posterior (blue). c) Absolute activating thresholds for the same 90 neurons in homogeneous space. d) Absolute activating thresholds for the same 90 neurons within the sample. The units in c,d) are given in terms of the base field strength of 100 V/m; one unit corresponds to 100 V/m. These results are given for the case of the unmyelinated axonal arbor and a negative biphasic activation pulse with the duration of 0.3 ms (Fig. S3). Download figure Open in new tab Fig. S6. a) Equivalent straight fiber density distribution of the sample characterizing its anisotropy. b) Activating threshold reduction factor for the first 90 neurons with the longest axonal arbor (sorted) for three polarizations of the primary electric field: medial-lateral (red), dorsal-ventral (green), and anterior-posterior (blue). c) Absolute activating thresholds for the same 90 neurons in homogeneous space. d) Absolute activating thresholds for the same 90 neurons within the sample. The units in c,d) are given in terms of the base field strength of 100 V/m; one unit corresponds to 100 V/m. These results are given for the case of the unmyelinated axonal arbor and a positive biphasic activation pulse with the duration of 0.3 ms (Fig. S3). Download figure Open in new tab Fig. S7. a) Equivalent straight fiber density distribution of the sample characterizing its anisotropy. b) Activating threshold reduction factor for the first 90 neurons with the longest axonal arbor (sorted) for three polarizations of the primary electric field: medial-lateral (red), dorsal-ventral (green), and anterior-posterior (blue). c) Absolute activating thresholds for the same 90 neurons in homogeneous space. d) Absolute activating thresholds for the same 90 neurons within the sample. The units in c,d) are given in terms of the base field strength of 100 V/m; one unit corresponds to 100 V/m. These results are given for the case of the myelinated axonal arbor and a negative rectangular activation monophasic pulse with the duration of 0.1 ms. Download figure Open in new tab Fig. S8. a) Equivalent straight fiber density distribution of the sample characterizing its anisotropy. b) Activating threshold reduction factor for the first 90 neurons with the longest axonal arbor (sorted) for three polarizations of the primary electric field: medial-lateral (red), dorsal-ventral (green), and anterior-posterior (blue). c) Absolute activating thresholds for the same 90 neurons in homogeneous space. d) Absolute activating thresholds for the same 90 neurons within the sample. The units in c,d) are given in terms of the base field strength of 100 V/m; one unit corresponds to 100 V/m. These results are given for the case of the myelinated axonal arbor and a positive rectangular activation monophasic pulse with the duration of 0.1 ms. Download figure Open in new tab Fig. S9. a) Equivalent straight fiber density distribution of the sample characterizing its anisotropy. b) Activating threshold reduction factor for the first 90 neurons with the longest axonal arbor (sorted) for three polarizations of the primary electric field: medial-lateral (red), dorsal-ventral (green), and anterior-posterior (blue). c) Absolute activating thresholds for the same 90 neurons in homogeneous space. d) Absolute activating thresholds for the same 90 neurons within the sample. The units in c,d) are given in terms of the base field strength of 100 V/m; one unit corresponds to 100 V/m. These results are given for the case of the myelinated axonal arbor and a negative biphasic activation pulse with the duration of 0.3 ms (Fig. S3). Download figure Open in new tab Fig. S10. a) Equivalent straight fiber density distribution of the sample characterizing its anisotropy. b) Activating threshold reduction factor for the first 90 neurons with the longest axonal arbor (sorted) for three polarizations of the primary electric field: medial-lateral (red), dorsal- ventral (green), and anterior-posterior (blue). c) Absolute activating thresholds for the same 90 neurons in homogeneous space. d) Absolute activating thresholds for the same 90 neurons within the sample. The units in c,d) are given in terms of the base field strength of 100 V/m; one unit corresponds to 100 V/m. These results are given for the case of the myelinated axonal arbor and a positive biphasic activation pulse with the duration of 0.3 ms (Fig. S3). Supplement D Estimates of absolute TMS activating threshold of mouse primary motor cortex from experimental data D1. Experimental procedure In vivo, TMS experiments were performed on male C57BL/6J mice (n¼6) using a custom rodent TMS coil with a magnetic core ( Fig. S11c ). Animals were anesthetized with sodium pentobarbital intraperitoneally (50 mg/kg), MEP (motor evoked potential) was recorded using a similar approach reported by Rotenberg et al. [ 13 ]. All procedures were approved by National Institute on Drug Abuse/NIH (NIDA) animal care and use committee. The coil was mounted to a customized three- axis micromanipulator. The focal electrical field point was carefully aligned to the targeted mouse motor cortex. Fine adjustments of the coil were made to induce limb twitch on the contralateral hindlimb, but not the ipsilateral hindlimb, any of the forelimbs or any other body part. More detailed description of TMS experiments is given in [ 14 ]. Detailed description of a similar custom focal TMS rodent coil was previously given in Ref. [ 15 ]. D2. Modeling setup Detailed brain model [ 16 ] ( Fig. S11a ) along with the model of the custom TMS rodent coil with a magnetic core ( Fig. S11c ) has been used to simulate the electric field distribution both in three principal planes ( Fig. S11b ) as well as just inside the brain compartments (Fig S11d). The skin shell conductivity was chosen as 0.3 S/m, the skull conductivity – 0.01 S/m. All brain compartments were assigned the conductivity of 0.3 S/m. Some conductivity variations were studied. They did not substantially alter the computational results. In general, the small mouse body quite significantly (by a factor of ∼2 or so) blocks the primary applied solenoidal electric field of the TMS coil. D3. Estimating electric field strength within mouse brain from measured motor evoked potentials Simulations were performed for two experimental parameters – the coil current of ∼3,000 A and dI/dt of 94 A/µs – corresponding to the induced MEPs. Fig. S11d shows the resulting magnitude of the total electric field just inside the brain compartments. The circle approximately labels the area of the primary motor cortex. One observes that the field strength there approaches 90 V/m. Conductivity variations of the head compartments lead to the field strength changes within the M1 area (white circle in Fig. S11 ) to within ±15% about this value. Similar results have been obtained for the total field magnitude distributions in three principal planes shown in Fig. S11e . Another mouse CAD model [ 17 ] has also been tested; it also generated very similar results. Download figure Open in new tab Fig. S11. a) CAD Mouse model [ 16 ]. b) Three principal planes for field output. c) Shape and positioning the custom TMS coil. d) Electric field magnitude distribution just inside mouse brain. Circle approximately indicates the primary motor cortex. D4. Comparison with previously reported computed data for mouse primary motor cortex In Ref. [ 17 ], the maximum fields strengths of 81.3 V/m,113.4 V/m, and 86.5 V/m in the brain have been computed for three different coil positions of a 25 mm figure-8 coil for the coil’s dI/dt of 100 A/μs In Ref. [ 18 ], it was concluded from numerical experiments that, for 100 A/μs, the predicted fields in mouse brain are on the order of 120 V/m. Our experimental data agree with these modeling predictions. They suggest that 90 V/m ± 15% is a reasonable MEP estimate in a mouse brain. Supplement E Verification of nested iterative algorithm E1. Test sample under study For comparison purposes, ten nearest neighbors of neuron #192 have been selected as shown in Fig. S12 below. The entire membrane mesh of the sample includes 16 million facets. After that, two methods have been applied: Global iterative BEM-FMM algorithm with the relative residual of 1e-4; Nested iterative BEM-FMM algorithm of the present study with the total number of 6 neuron-to-neuron outer iterations (as shown in Fig. 3 of the main text) and with the relative residual of 1e-4 for the inner iterative solution. The y-polarized primary electric field of 100 V/m (corresponding to a vertical direction in Fig. S12 ) has been applied in both cases. The global algorithm converges after 15 iterations, in ∼1,000 seconds, when using a 3 GHz workstation. Download figure Open in new tab Fig. S12. Truncated sample mesh of 11 neighbors of neuron #192. E2. Comparison results The two solutions have been compared to each other with regard to the induced membrane charge densities. Fig. S13 shows the relative 2-norm error in the local charge density for every neuron (11 in total) separately. One observes that this error does not exceed 1% on average. Download figure Open in new tab Fig. S13. Relative 2-norm error for the induced charge density between the two solutions for every individual neuron of the sample. Supplement F Computer code for cathodic stimulation F1. Computer code for cathodic stimulation + noise with arbitrarily varying parameters This code supports Section 4.3 (Discussion Section) of the main text. It is written within MATLAB platform and should run on Windows machines. A simple numerical experiment with a Hodgkin-Huxley axon originally studied by Rattay [ 19 ] is considered. A cathodic stimulation of the 10 mm long axon by a point electrode located at the origin is modeled. The default axonal radius is 5 µm, the length constant λ is 1.5 mm; the time constant τ is 3.3 ms, the resting potential is – 65 mV. A monophasic electric pulse (T=0.1 ms) is applied. Field perturbers are 2D cylinders with non-conducting walls located perpendicular to the axon and with a large radius of 200 µm, chosen for better visualization. An analytical solution for every infinite 2D cylinder secondary field is employed (cf., for example, [ 20 ]). First, the MATLAB script “script01_define_axon_HH_model.m” should be executed. There, the major biophysical and computational parameters are defined. All units are SI units. Extensive comments are provided for every quantity of interest. Second, the MATLAB script “script02_obtain_axon_solution.m” should be executed. This script runs biophysical computations for the transmembrane potential and the ionic transmembrane current density. If necessary, the total transmembrane current density could be plotted as well. Two biophysical solutions are executed simultaneously: the first one is for the homogeneous medium (red color) and the second one is for the perturbed extracellular field (blue color). Fig. 9 of the main text shows the graphical output of the MATLAB script “script02_obtain_axon_solution.m” that is updated at every time moment. Solution parameters (electrode current, axon data, time course, spatial and temporal resolutions, cylinder radii, their spacing and positions) as well as the output format can be modified using both executable MATLAB scripts. Other MATLAB functions should not be changed. F2. Code availability The Dropbox link [ 21 ] contains a fully executable computational module built within the MATLAB platform under Windows and a short guide that replicates the above example. This enables both testing of different polarization directions as well as extraction of collinear fields along any desired axonal path. Acknowledgements Authors express their gratitude to Dr. Bethanny Danskin and the Virtual Observatory of the Cortex (VORTEX) program of the Allen Institute for Brain Science for their assistance in processing the L2/L3 brain sample, to Dr. Clayton S. Bingham of NINDS/NIH and Dr. Cameron C. McIntyre of Duke University for shaping early stages of this study, to Dr. Ted Carnevale of Yale School of Medicine for his guidance in utilizing NEURON software with superior spatial resolution, and to Dr. Padma Sundaram of Massachusetts General Hospital for useful discussions. Special thanks are extended to Dr. Ermal Toto and Dr. James Kingsley of Academic and Research Computing at Worcester Polytechnic Institute for their great support in assembling and executing highly demanding parallel computational tasks. This study has received support from the NIMH grant R01MH130490 and NIBIB grant R01EB035484 (SNM), NIDCD grant R01DC020891 (ARN), NINDS grant U24NS120053 (BD) as well as the Intramural Research Program of NIH, NIMH ZIAMH002955 (Z-DD), and Intramural Research Program of NIH, NIDA ZIADA000638 (HL). Finally, we express our sincerest thanks to the two anonymous reviewers for their dedication and comprehensive, professional comments, which shaped both the methods and conclusions of this study. We are also thankful to the deputy editor for their effort. Footnotes Significant changes to the text, including revisions to the methods, results, and discussion. New figures have been added to clarify the methods and results. Important limitations to the study have been identified. References [1]. ↵ Siebner HR , Funke K , Aberra AS , Antal A , Bestmann S , Chen R , Classen J , Davare M , Di Lazzaro V , Fox PT , Hallett M , Karabanov AN , Kesselheim J , Beck MM , Koch G , Liebetanz D , Meunier S , Miniussi C , Paulus W , Peterchev AV , Popa T , Ridding MC , Thielscher A , Ziemann U , Rothwell JC , Ugawa Y . Transcranial magnetic stimulation of the brain: What is stimulated? - A consensus and critical position paper . Clin Neurophysiol . 2022 Aug ; 140 : 59 – 97 . doi: 10.1016/j.clinph.2022.04.022 . OpenUrl CrossRef PubMed [2]. ↵ Liu A , Vöröslakos M , Kronberg G , Henin S , Krause MR , Huang Y , Opitz A , Mehta A , Pack CC , Krekelberg B , Berényi A , Parra LC , Melloni L , Devinsky O , Buzsáki G . Immediate neurophysiological effects of transcranial electrical stimulation . Nat Commun . 2018 Nov 30; 9 ( 1 ): 5092 . doi: 10.1038/s41467-018-07233-7 . OpenUrl CrossRef PubMed [3]. ↵ Lozano AM , Lipsman N , Bergman H , Brown P , Chabardes S , Chang JW , Matthews K , McIntyre CC , Schlaepfer TE , Schulder M , Temel Y , Volkmann J , Krauss JK . Deep brain stimulation: current challenges and future directions . Nat Rev Neurol . 2019 Mar ; 15 ( 3 ): 148 – 160 . doi: 10.1038/s41582-018-0128-2 . OpenUrl CrossRef PubMed [4]. ↵ Deng ZD , Argyelan M , Miller J , Quinn DK , Lloyd M , Jones TR , Upston J , Erhardt E , McClintock SM , Abbott CC . Electroconvulsive therapy, electric field, neuroplasticity, and clinical outcomes . Mol Psychiatry . 2022 Mar ; 27 ( 3 ): 1676 – 1682 . doi: 10.1038/s41380-021-01380-y . OpenUrl CrossRef PubMed [5]. ↵ Turner NL , Macrina T , Bae JA , Yang R , Wilson AM , Schneider-Mizell C , Lee K , Lu R , Wu J , Bodor AL , Bleckert AA , Brittain D , Froudarakis E , Dorkenwald S , Collman F , Kemnitz N , Ih D , Silversmith WM , Zung J , Zlateski A , Tartavull I , Yu SC , Popovych S , Mu S , Wong W , Jordan CS , Castro M , Buchanan J , Bumbarger DJ , Takeno M , Torres R , Mahalingam G , Elabbady L , Li Y , Cobos E , Zhou P , Suckow S , Becker L , Paninski L , Polleux F , Reimer J , Tolias AS , Reid RC, da Costa NM , Seung HS. Reconstruction of neocortex: Organelles, compartments, cells, circuits, and activity. Cell . 2022 Mar 17; 185 ( 6 ): 1082 – 1100 .e24. doi: 10.1016/j.cell.2022.01.023 . OpenUrl CrossRef [6]. ↵ Dorkenwald S , Turner NL , Macrina T , Lee K , Lu R , Wu J , Bodor AL , Bleckert AA , Brittain D , Kemnitz N , Silversmith WM , Ih D , Zung J , Zlateski A , Tartavull I , Yu SC , Popovych S , Wong W , Castro M , Jordan CS , Wilson AM , Froudarakis E , Buchanan J , Takeno MM , Torres R , Mahalingam G , Collman F , Schneider-Mizell CM , Bumbarger DJ , Li Y , Becker L , Suckow S , Reimer J , Tolias AS , Macarico da Costa N , Reid RC , Seung HS . Binary and analog variation of synapses between cortical pyramidal neurons . Elife . 2022 Nov 16; 11 : e76120 . doi: 10.7554/eLife.76120 . OpenUrl CrossRef PubMed [7]. ↵ MICrONs Consortium et al. Functional connectomics spanning multiple areas of mouse visual cortex . bioRxiv originally published 2021.07.28.454025/updated 2023.04.19; doi: https://www.biorxiv.org/content/10.1101/2021.07.28.454025v3 . [8]. ↵ Eisenstein M . A milestone map of mouse-brain connectivity reveals challenging new terrain for scientists . Nature . 2024 Apr ; 628 (8008): 677 -679. doi: 10.1038/d41586-024-01096-3 . OpenUrl CrossRef PubMed [9]. ↵ Shapson-Coe A , Januszewski M , Berger DR , Pope A , Wu Y , Blakely T , Schalek RL , Li PH , Wang S , Maitin-Shepard J , Karlupia N , Dorkenwald S , Sjostedt E , Leavitt L , Lee D , Troidl J , Collman F , Bailey L , Fitzmaurice A , Kar R , Field B , Wu H , Wagner-Carena J , Aley D , Lau J , Lin Z , Wei D , Pfister H , Peleg A , Jain V , Lichtman JW . A petavoxel fragment of human cerebral cortex reconstructed at nanoscale resolution . Science . 2024 May 10; 384 ( 6696 ):eadk4858. doi: 10.1126/science.adk4858 . OpenUrl CrossRef [10]. ↵ Aberra AS , Peterchev AV , Grill WM . Biophysically realistic neuron models for simulation of cortical stimulation . J Neural Eng . 2018 Dec ; 15 ( 6 ): 066023 . doi: 10.1088/1741-2552/aadbb1 . OpenUrl CrossRef PubMed [11]. ↵ Aberra AS , Wang B , Grill WM , Peterchev AV . Simulation of transcranial magnetic stimulation in head model with morphologically-realistic cortical neurons . Brain Stimul . 2020 Jan -Feb; 13 ( 1 ): 175 – 189 . doi: 10.1016/j.brs.2019.10.002 . OpenUrl CrossRef PubMed [12]. ↵ Weise K , Worbs T , Kalloch B , Souza VH , Jaquier AT , Van Geit W , Thielscher A , Knösche TR . Directional Sensitivity of Cortical Neurons Towards TMS Induced Electric Fields . Imaging Neuroscience ( 2023 ) 1 : 1 – 22 . doi: 10.1162/imag_a_00036 . OpenUrl CrossRef [13]. ↵ Bingham CS , McIntyre CC . Subthalamic deep brain stimulation of an anatomically detailed model of the human hyperdirect pathway . J Neurophysiol . 2022 May 1; 127 ( 5 ): 1209 – 1220 . doi: 10.1152/jn.00004.2022 . OpenUrl CrossRef PubMed [14]. ↵ Noor MS , McIntyre CC . Biophysical characterization of local field potential recordings from directional deep brain stimulation electrodes . Clin Neurophysiol . 2021 Jun ; 132 ( 6 ): 1321 – 1329 . doi: 10.1016/j.clinph.2021.01.027 . OpenUrl CrossRef [15]. ↵ Bower KL , McIntyre CC . Deep brain stimulation of terminating axons . Brain Stimul . 2020 Nov Dec;1 3 ( 6 ): 1863 – 1870 . doi: 10.1016/j.brs.2020.09.001 . OpenUrl CrossRef 16. ↵ McIntyre CC . Private communication . Oct. 2022 . [17]. ↵ Liu , W.K. , Li , S. & Park , H.S . Eighty Years of the Finite Element Method: Birth , Evolution, and Future. Arch Computat Methods Eng 29 , 4431 – 4453 ( 2022 ). doi: 10.1007/s11831-022-09740-9 . OpenUrl CrossRef [18]. ↵ Ye H , Steiger A . Neuron matters: electric activation of neuronal tissue is dependent on the interaction between the neuron and the electric field . J Neuroeng Rehabil . 2015 Aug 12; 12 : 65 . doi: 10.1186/s12984-015-0061-1 . OpenUrl CrossRef PubMed [19]. ↵ Maccabee PJ , Amassian VE , Eberle LP , Cracco RQ . Magnetic coil stimulation of straight and bent amphibian and mammalian peripheral nerve in vitro: locus of excitation . J Physiol . 1993 Jan ; 460 : 201 – 19 . doi: 10.1113/jphysiol.1993.sp019467 . PMID: 8487192 ; PMCID: PMC1175209 . OpenUrl CrossRef PubMed Web of Science [20]. ↵ Joucla S , Branchereau P , Cattaert D , Yvert B . Extracellular neural microstimulation may activate much larger regions than expected by simulations: a combined experimental and modeling study . PLoS One . 2012 ; 7 ( 8 ): e41324 . doi: 10.1371/journal.pone.0041324 . OpenUrl CrossRef PubMed [21]. ↵ Fellner A , Heshmat A , Werginz P , Rattay F . A finite element method framework to model extracellular neural stimulation . J Neural Eng . 2022 Apr 7; 19 ( 2 ). doi: 10.1088/1741-2552/ac6060 . OpenUrl CrossRef [22]. ↵ Agudelo-Toro A , Neef A . Computationally efficient simulation of electrical activity at cell membranes interacting with self-generated and externally imposed electric fields . J Neural Eng . 2013 Apr ; 10 ( 2 ): 026019 . doi: 10.1088/1741-2560/10/2/026019 . OpenUrl CrossRef PubMed [23]. ↵ Meffin H , Tahayori B , Sergeev EN , Mareels IM , Grayden DB , Burkitt AN . Modelling extracellular electrical stimulation: part 3. Derivation and interpretation of neural tissue equations . J Neural Eng . 2014 Dec ; 11 ( 6 ): 065004 . doi: 10.1088/1741-2560/11/6/065004 . OpenUrl CrossRef PubMed [24]. ↵ Xylouris K , Wittum G . A three-dimensional mathematical model for the signal propagation on a neuron’s membrane . Front Comput Neurosci . 2015 Jul 17; 9 : 94 . doi: 10.3389/fncom.2015.00094 . OpenUrl CrossRef [25]. ↵ Monfared O , Tahayori B , Freestone D , Nešić D , Grayden DB , Meffin H . Determination of the electrical impedance of neural tissue from its microscopic cellular constituents . J Neural Eng . 2020 Jan 24; 17 ( 1 ): 016037 . doi: 10.1088/1741-2552/ab560a . OpenUrl CrossRef PubMed [26]. ↵ Pelot NA , Behrend CE , Grill WM . On the parameters used in finite element modeling of compound peripheral nerves . J Neural Eng . 2019 Feb ; 16 ( 1 ): 016007 . doi: 10.1088/1741-2552/aaeb0c . OpenUrl CrossRef PubMed [27]. ↵ Greengard L , Rokhlin V . A fast algorithm for particle simulations . J. Comput. Phys . 1987 ; 73 ( 2 ): 325 – 348 . doi: 10.1016/0021-9991(87)90140-9 . OpenUrl CrossRef Web of Science [28]. ↵ Greengard L. , Rokhlin V . A new version of the Fast Multipole Method for the Laplace equation in three dimensions . Acta Numerica , 6 : 229 – 269 , 1997 . OpenUrl CrossRef [29]. ↵ Gimbutas Z , Greengard L , Magland J , Rachh M , Rokhlin V . fmm3D Documentation . Release 0.1.0 . 2019– 2024 . Online: https://github.com/flatironinstitute/FMM3D & https://github.com/flatironinstitute/FMM3D/blob/master/fmm3d_manual.pdf . [30]. ↵ Greengard L , Gueyffier D , Martinsson PG , Rokhlin V ( May 2009 ). “ Fast direct solvers for integral equations in complex three-dimensional domains ”. Acta Numerica . 18 : 243 – 275 . doi: 10.1017/S0962492906410011 . ISSN 1474-0508. S2CID 58895952. OpenUrl CrossRef [31]. ↵ Makarov SN , Noetscher GM , Raij T , Nummenmaa A . A Quasi-Static Boundary Element Approach with Fast Multipole Acceleration for High-Resolution Bioelectromagnetic Models . IEEE Trans Biomed Eng . 2018 Mar 7. doi: 10.1109/TBME.2018.2813261 . OpenUrl CrossRef PubMed [32]. ↵ Noetscher GM , Tang D , Nummenmaa AR , Bingham CS , McIntyre CC , Makaroff SN . Estimations of Charge Deposition onto Convoluted Axon Surfaces within Extracellular Electric Fields . IEEE Trans Biomed Eng . 2023 Aug 3; doi: 10.1109/TBME.2023.3299734 . OpenUrl CrossRef [33]. ↵ Makaroff SN , Nummenmaa AR , Noetscher GM , Qi Z , McIntyre CC , Bingham CS . Influence of charges deposited on membranes of human hyperdirect pathway axons on depolarization during subthalamic deep brain stimulation . J Neural Eng . 2023 Jul 19; 20 ( 4 ). doi: 10.1088/1741-2552/ace5de . OpenUrl CrossRef [34]. ↵ Krassowska W , Neu JC . Response of a single cell to an external electric field . Biophys J . 1994 Jun ; 66 ( 6 ): 1768 – 76 . doi: 10.1016/S0006-3495(94)80971-3 . OpenUrl CrossRef PubMed Web of Science [35]. Khadka N , Bikson M. Neurocapillary-Modulation . Neuromodulation . 2022 Dec ; 25 ( 8 ): 1299 – 1311 . doi: 10.1111/ner.13338 . OpenUrl CrossRef PubMed [36]. Khadka N , Poon C , Cancel LM , Tarbell JM , Bikson M . Multi-scale multi-physics model of brain interstitial water flux by transcranial Direct Current Stimulation . J Neural Eng . 2023 Jul 24; 20 ( 4 ):10.1088/1741-2552/ace4f4. doi: 10.1088/1741-2552/ace4f4 . OpenUrl CrossRef [37]. ↵ Kress R 1999 Linear Integral Equations 2nd ed, New York , Springer . [38]. ↵ Ponasso GN . A survey on integral equations for bioelectric modeling . Phys Med Biol . 2024 Aug 28; 69 ( 17 ):10.1088/1361-6560/ad66a9. doi: 10.1088/1361-6560/ad66a9 . OpenUrl CrossRef [39]. ↵ Onorato I , D’Alessandro G , Di Castro MA , Renzi M , Dobrowolny G , Musarò A , Salvetti M , Limatola C , Crisanti A , Grassi F . Noise Enhances Action Potential Generation in Mouse Sensory Neurons via Stochastic Resonance . PLoS One . 2016 Aug 15; 11 ( 8 ): e0160950 . doi: 10.1371/journal.pone.0160950 . OpenUrl CrossRef PubMed [40]. ↵ van der Groen O , Potok W , Wenderoth N , Edwards G , Mattingley JB , Edwards D . Using noise for the better: The effects of transcranial random noise stimulation on the brain and behavior . Neurosci Biobehav Rev . 2022 Jul ; 138 : 104702 . doi: 10.1016/j.neubiorev.2022.104702 . OpenUrl CrossRef PubMed [41]. ↵ Rattay F . Analysis of models for external stimulation of axons . IEEE Trans Biomed Eng . 1986 Oct ; 33 ( 10 ): 974 – 7 . doi: 10.1109/TBME.1986.325670 . PMID: 3770787 . OpenUrl CrossRef PubMed Web of Science [42]. ↵ Rattay F . Modeling the excitation of fibers under surface electrodes . IEEE Trans Biomed Eng . 1988 Mar ; 35 ( 3 ): 199 – 202 . doi: 10.1109/10.1362 . OpenUrl CrossRef PubMed [43]. ↵ Rattay F . Analysis of models for extracellular fiber stimulation . IEEE Trans Biomed Eng . 1989 Jul ; 36 ( 7 ): 676 – 82 . doi: 10.1109/10.32099 . OpenUrl CrossRef PubMed Web of Science [44]. ↵ Czerwonky DM , Aberra AS , Gomez LJ . A Boundary Element Method of Bidomain Modeling for Predicting Cellular Responses to Electromagnetic Fields . J Neural Eng . 2024 Jun 25; 21 ( 3 ). doi: 10.1088/1741-2552/ad5704 . OpenUrl CrossRef [45]. ↵ Hines ML , Carnevale NT . The NEURON simulation environment . Neural Comput . 1997 Aug 15; 9 ( 6 ): 1179 – 209 . doi: 10.1162/neco.1997.9.6.1179 . PMID: 9248061 . OpenUrl CrossRef PubMed Web of Science [46]. ↵ Markram H , Muller E , Ramaswamy S , Reimann MW , Abdellah M , Sanchez CA , Ailamaki A , Alonso-Nanclares L , Antille N , Arsever S , Kahou GA , Berger TK , Bilgili A , Buncic N , Chalimourda A , Chindemi G , Courcol JD , Delalondre F , Delattre V , Druckmann S , Dumusc R , Dynes J , Eilemann S , Gal E , Gevaert ME , Ghobril JP , Gidon A , Graham JW , Gupta A , Haenel V , Hay E , Heinis T , Hernando JB , Hines M , Kanari L , Keller D , Kenyon J , Khazen G , Kim Y , King JG , Kisvarday Z , Kumbhar P , Lasserre S , Le Bé JV , Magalhães BR , Merchán-Pérez A , Meystre J , Morrice BR , Muller J , Muñoz-Céspedes A , Muralidhar S , Muthurasa K , Nachbaur D , Newton TH , Nolte M , Ovcharenko A , Palacios J , Pastor L , Perin R , Ranjan R , Riachi I , Rodríguez JR , Riquelme JL , Rössert C , Sfyrakis K , Shi Y , Shillcock JC , Silberberg G , Silva R , Tauheed F , Telefont M , Toledo-Rodriguez M , Tränkler T , Van Geit W , Díaz JV , Walker R , Wang Y , Zaninetta SM , DeFelipe J , Hill SL , Segev I , Schürmann F . Reconstruction and Simulation of Neocortical Microcircuitry . Cell . 2015 Oct 8; 163 ( 2 ): 456 – 92 . doi: 10.1016/j.cell.2015.09.029 . OpenUrl CrossRef PubMed [47]. ↵ Makarov SN , Noetscher GM , Nazarian A . Low-Frequency Electromagnetic Modeling of Electrical and Biological Systems Using MATLAB . New York : Wiley , 2015 , 598 p., ISBN: 978-1119052562. [48]. ↵ Pilgrim C , Reisert I , Grab D . Volume densities and specific surfaces of neuronal and glial tissue elements in the rat supraoptic nucleus . J Comp Neurol . 1982 Nov 10; 211 ( 4 ): 427 – 31 . doi: 10.1002/cne.902110409 . OpenUrl CrossRef PubMed [49]. ↵ Lehmenkühler A , Syková E , Svoboda J , Zilles K , Nicholson C . Extracellular space parameters in the rat neocortex and subcortical white matter during postnatal development determined by diffusion analysis . Neuroscience . 1993 Jul ; 55 ( 2 ): 339 – 51 . doi: 10.1016/0306-4522(93)90503-8 . OpenUrl CrossRef PubMed Web of Science [50]. ↵ Mota B , Herculano-Houzel S . All brains are made of this: a fundamental building block of brain matter with matching neuronal and glial masses . Front Neuroanat . 2014 Nov 12; 8 : 127 . doi: 10.3389/fnana.2014.00127 . OpenUrl CrossRef [51]. Clark J , Plonsey R . A mathematical evaluation of the core conductor model . Biophys J . 1966 Jan ; 6 ( 1 ): 95 – 112 . doi: 10.1016/S0006-3495(66)86642-0 . OpenUrl CrossRef PubMed Web of Science [52]. Clark J , Plonsey R . The extracellular potential field of the single active nerve fiber in a volume conductor . Biophys J . 1968 Jul ; 8 ( 7 ): 842 – 64 . doi: 10.1016/S0006-3495(68)86524-5 . OpenUrl CrossRef PubMed Web of Science [53]. Klee M , Plonsey R . Stimulation of spheroidal cells--the role of cell shape . IEEE Trans Biomed Eng . 1976 Jul ; 23 ( 4 ): 347 – 54 . OpenUrl PubMed Web of Science [54]. Huang Y , Liu AA , Lafon B , Friedman D , Dayan M , Wang X , Bikson M , Doyle WK , Devinsky O , Parra LC . Measurements and models of electric fields in the in vivo human brain during transcranial electric stimulation . Elife . 2017 Feb 7; 6 : e18834 . doi: 10.7554/eLife.18834 OpenUrl CrossRef PubMed [55]. ↵ Pourtaheri N , Ying W , Kim JM , Henriquez CS . Thresholds for transverse stimulation: fiber bundles in a uniform field . IEEE Trans Neural Syst Rehabil Eng . 2009 Oct ; 17 ( 5 ): 478 – 86 . doi: 10.1109/TNSRE.2009.2033424 . OpenUrl CrossRef PubMed [56]. ↵ Klee M , Plonsey R . Stimulation of spheroidal cells-The role of cell shape . 1976 . IEEE Trans Biomed Eng . Jul; 23 ( 4 ): 347 – 54 . PMID: 1278928 . OpenUrl PubMed Web of Science [57]. ↵ Pavlin M , Kanduser M , Rebersek M , Pucihar G , Hart FX , Magjarevic R , Miklavcic D . Effect of cell electroporation on the conductivity of a cell suspension . Biophys J . 2005 Jun ; 88 ( 6 ): 4378 – 90 . doi: 10.1529/biophysj.104.048975 . OpenUrl CrossRef PubMed [58]. ↵ Pavlin M , Slivnik T , Miklavcic D . Effective conductivity of cell suspensions . IEEE Trans Biomed Eng . 2002 Jan ; 49 ( 1 ): 77 – 80 . doi: 10.1109/10.972843 . OpenUrl CrossRef PubMed [59]. ↵ Maxwell , J. C . 1873 . Treatise on Electricity and Magnetism . Oxford University Press , London, UK . [60]. ↵ Saad Y . Iterative Methods for Sparse Linear Systems . 2nd Ed ., Society for Industrial and Applied Mathematics . 2003 . ISBN 978-0-89871-534-7. [61]. ↵ Siebner HR , Funke K , Aberra AS , Antal A , Bestmann S , Chen R , Classen J , Davare M , Di Lazzaro V , Fox PT , Hallett M , Karabanov AN , Kesselheim J , Beck MM , Koch G , Liebetanz D , Meunier S , Miniussi C , Paulus W , Peterchev AV , Popa T , Ridding MC , Thielscher A , Ziemann U , Rothwell JC , Ugawa Y . Transcranial magnetic stimulation of the brain: What is stimulated? - A consensus and critical position paper . Clin Neurophysiol . 2022 Aug ; 140 : 59 – 97 . doi: 10.1016/j.clinph.2022.04.022 . OpenUrl CrossRef PubMed [62]. MICrONS Consortium . Online: https://www.microns-explorer.org/cortical-mm3 [63]. Jia Wan , Wanhua Li , Atmadeep Banerjee , Jason Ken Adhinarta , Evelina Sjostedt , Jingpeng Wu , Jeff Lichtman , Hanspeter Pfister , Donglai Wei . TriSAM: Tri-Plane SAM for zero-shot cortical blood vessel segmentation in VEM images . arXiv :2401.13961v1 . Jan. 2024 . doi: 10.48550/arXiv.2401.13961 . OpenUrl CrossRef [64]. ↵ Zhao Z , Shirinpour S , Tran H , Wischnewski M , Opitz A . Intensity- and frequency-specific effects of transcranial alternating current stimulation are explained by network dynamics . J Neural Eng . 2024 Apr 03; 21 ( 2 ). doi: 10.1088/1741-2552/ad37d9 . OpenUrl CrossRef [65]. Makaroff et al. 2024 [Data set]. Electromagnetic modeling data for the IARPA MICrONS Pinky100 mouse brain dataset . BossDB Archive . doi https://10.60533/boss-2024-vf8f . Online: https://bossdb.org/project/makaroff2024 [66]. ↵ Bingham CS , McIntyre CC . Subthalamic deep brain stimulation of an anatomically detailed model of the human hyperdirect pathway . J Neurophysiol , vol. 127 , no. 5 , pp. 1209 – 1220 , May 2022 . OpenUrl CrossRef PubMed [67]. ↵ Noor MS , McIntyre CC . Biophysical characterization of local field potential recordings from directional deep brain stimulation electrodes . Clin Neurophysiol , vol. 132 , no. 6 , pp. 13211329, Jun. 2021 .. [68]. ↵ Bower KL , McIntyre CC . Deep brain stimulation of terminating axons . Brain Stimul , vol. 13 , no. 6 , pp. 1863 – 1870 , Nov. 2020 . OpenUrl CrossRef PubMed [69]. ↵ Musk E ; Neuralink. An Integrated Brain-Machine Interface Platform With Thousands of Channels . J Med Internet Res . 2019 Oct 31; 21 ( 10 ): e16194 . doi: 10.2196/16194 . PMID: 31642810 ; PMCID: PMC6914248 . OpenUrl CrossRef PubMed [70]. ↵ Saha S , Mamun KA , Ahmed K , Mostafa R , Naik GR , Darvishi S , Khandoker AH , Baumert M . Progress in Brain Computer Interface: Challenges and Opportunities . Front Syst Neurosci . 2021 Feb 25; 15 : 578875 . doi: 10.3389/fnsys.2021.578875 . OpenUrl CrossRef [71]. ↵ Card NS , Wairagkar M , Iacobacci C , Hou X , Singer-Clark T , Willett FR , Kunz EM , Fan C , Vahdati Nia M , Deo DR , Srinivasan A , Choi EY , Glasser MF , Hochberg LR , Henderson JM , Shahlaie K , Stavisky SD , Brandman DM . An Accurate and Rapidly Calibrating Speech Neuroprosthesis . N Engl J Med . 2024 Aug 15; 391 ( 7 ): 609 – 618 . doi: 10.1056/NEJMoa2314132 . OpenUrl CrossRef PubMed [72]. ↵ Computer code to Supplement F . 2024 . Dropbox Online . References – Supplements A, B, C, D, E, F [1]. ↵ [ MICrONs Consortium et al. Functional connectomics spanning multiple areas of mouse visual cortex . bioRxiv originally published 2021.07.28.454025/updated 2023.04.19; doi: https://www.biorxiv.org/content/10.1101/2021.07.28.454025v3 . Online: https://www.microns-explorer.org/phase1 . Release Notes (v185). [2]. ↵ CGAL online: https://github.com/CGAL/cgal [3]. ↵ Pierre Alliez and David Cohen-Steiner and Michael Hemmer and Cédric Portaneri and Mael Rouxel-Labbé . 3D Alpha Wrapping. In CGAL User and Reference Manual . CGAL Editorial Board, 5.6 edition , 2023 . [4]. ↵ Pilgrim C , Reisert I , Grab D Volume densities and specific surfaces of neuronal and glial tissue elements in the rat supraoptic nucleus J Comp Neurol 1982 Nov 10; 211 ( 4 ): 427 – 31 . doi: 10.1002/cne.902110409 . OpenUrl CrossRef PubMed [5]. ↵ Lehmenkühler A , Syková E , Svoboda J , Zilles K , Nicholson C Extracellular space parameters in the rat neocortex and subcortical white matter during postnatal development determined by diffusion analysis Neuroscience 1993 Jul ; 55 ( 2 ): 339 – 51 . doi: 10.1016/0306-4522(93)90503-8 . OpenUrl CrossRef PubMed Web of Science [6]. ↵ Mota B , Herculano-Houzel S All brains are made of this: a fundamental building block of brain matter with matching neuronal and glial masses Front Neuroanat 2014 Nov 12; 8 : 127 . doi: 10.3389/fnana.2014.00127 . OpenUrl CrossRef [7]. ↵ Keller D , Erö C , Markram H Cell Densities in the Mouse Brain: A Systematic Review Front Neuroanat 2018 Oct 23; 12 : 83 . doi: 10.3389/fnana.2018.00083 . OpenUrl CrossRef PubMed [8]. ↵ Sven Holcombe ( 2024 ). inpolyhedron - are points inside a triangulated volume? ( https://www.mathworks.com/matlabcentral/fileexchange/37856-inpolyhedron-are-points-inside-a-triangulated-volume ), MATLAB Central File Exchange. Retrieved February 5, 2024. [9]. ↵ Hines ML , Carnevale NT The NEURON simulation environment Neural Comput 1997 Aug 15; 9 ( 6 ): 1179 – 209 . doi: 10.1162/neco.1997.9.6.1179 . PMID: 9248061 . OpenUrl CrossRef PubMed Web of Science [10]. ↵ Aberra AS , Peterchev AV , Grill WM Biophysically realistic neuron models for simulation of cortical stimulation J Neural Eng 2018 Dec ; 15 ( 6 ): 066023 . doi: 10.1088/1741-2552/aadbb1 . OpenUrl CrossRef PubMed [11]. ↵ Aberra AS , Wang B , Grill WM , Peterchev AV Simulation of transcranial magnetic stimulation in head model with morphologically-realistic cortical neurons Brain Stimul 2020 Jan -Feb; 13 ( 1 ): 175 – 189 . doi: 10.1016/j.brs.2019.10.002 . OpenUrl CrossRef PubMed [12]. ↵ Weise K , Worbs T , Kalloch B , Souza VH , Jaquier AT , Van Geit W , Thielscher A , Knösche TR Directional Sensitivity of Cortical Neurons Towards TMS Induced Electric Fields bioRxiv 2023.07.06.547913; doi: 10.1101/2023.07.06.547913 . OpenUrl Abstract / FREE Full Text [13]. ↵ Rotenberg A , Muller PA , Vahabzadeh-Hagh AM , Navarro X , López-Vales R , Pascual- Leone A , Jensen F Lateralization of forelimb motor evoked potentials by transcranial magnetic stimulation in rats Clin Neurophysiol 2010 Jan ; 121 ( 1 ): 104 – 8 . doi: 10.1016/j.clinph.2009.09.008 . Epub 2009 Nov 8. PMID: 19900839 ; PMCID: PMC2818443 . OpenUrl CrossRef PubMed Web of Science [14]. ↵ Meng Q , Jing L , Badjo JP , Du X , Hong E , Yang Y , Lu H , Choa FS A novel transcranial magnetic stimulator for focal stimulation of rodent brain Brain Stimul 2018 May - Jun; 11 ( 3 ): 663 – 665 . doi: 10.1016/j.brs.2018.02.018 . OpenUrl CrossRef PubMed [15]. ↵ Makaroff SN , Nguyen H , Meng Q , Lu H , Nummenmaa AR , Deng ZD Modeling transcranial magnetic stimulation coil with magnetic cores J Neural Eng 2023 Jan 25; 20 ( 1 ): 016028 . doi: 10.1088/1741-2552/acae0d . PMID: 36548994 ; PMCID: PMC10481791 . OpenUrl CrossRef PubMed [16]. ↵ Serano P , Makaroff S , Ackerman JL , Nummenmaa A , Noetscher GM Detailed high quality surface-based mouse CAD model suitable for electromagnetic simulations Biomed Phys Eng Express 2023 Nov 30; 10 ( 1 ):10.1088/2057-1976/ad0e14. doi: 10.1088/2057-1976/ad0e14 . OpenUrl CrossRef [17]. ↵ Alekseichuk I , Mantell K , Shirinpour S , Opitz A Comparative modeling of transcranial magnetic and electric stimulation in mouse, monkey, and human Neuroimage 2019 Jul 1; 194 : 136 – 148 . doi: 10.1016/j.neuroimage.2019.03.044 . OpenUrl CrossRef PubMed [18]. ↵ Koponen LM , Stenroos M , Nieminen JO , Jokivarsi K , Gröhn O , Ilmoniemi RJ Individual head models for estimating the TMS-induced electric field in rat brain Sci Rep 2020 Oct 15; 10 ( 1 ): 17397 . doi: 10.1038/s41598-020-74431-z . OpenUrl CrossRef PubMed [19]. ↵ Rattay F Analysis of models for extracellular fiber stimulation IEEE Trans Biomed Eng 1989 Jul ; 36 ( 7 ): 676 – 82 . doi: 10.1109/10.32099 . OpenUrl CrossRef PubMed Web of Science [20]. ↵ Noetscher GM , Tang D , Nummenmaa AR , Bingham CS , McIntyre CC , Makaroff SN Estimations of Charge Deposition onto Convoluted Axon Surfaces within Extracellular Electric Fields IEEE Trans Biomed Eng 2023 Aug 3; doi: 10.1109/TBME.2023.3299734 . OpenUrl CrossRef [21]. ↵ Computer codes to Supplement F 2024 Dropbox Online . View the discussion thread. Back to top Previous Next Posted November 30, 2024. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Enabling Electric Field Model of Microscopically Realistic Brain Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Enabling Electric Field Model of Microscopically Realistic Brain Zhen Qi , Gregory M. Noetscher , Alton Miles , Konstantin Weise , Thomas R. Knösche , Cameron R. Cadman , Alina R. Potashinsky , Kelu Liu , William A. Wartman , Guillermo Nunez Ponasso , Marom Bikson , Hanbing Lu , Zhi-De Deng , Aapo R. Nummenmaa , Sergey N. Makaroff bioRxiv 2024.04.04.588004; doi: https://doi.org/10.1101/2024.04.04.588004 Share This Article: Copy Citation Tools Enabling Electric Field Model of Microscopically Realistic Brain Zhen Qi , Gregory M. Noetscher , Alton Miles , Konstantin Weise , Thomas R. Knösche , Cameron R. Cadman , Alina R. Potashinsky , Kelu Liu , William A. Wartman , Guillermo Nunez Ponasso , Marom Bikson , Hanbing Lu , Zhi-De Deng , Aapo R. Nummenmaa , Sergey N. Makaroff bioRxiv 2024.04.04.588004; doi: https://doi.org/10.1101/2024.04.04.588004 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Neuroscience Subject Areas All Articles Animal Behavior and Cognition (7644) Biochemistry (17728) Bioengineering (13916) Bioinformatics (42037) Biophysics (21488) Cancer Biology (18636) Cell Biology (25552) Clinical Trials (138) Developmental Biology (13401) Ecology (19940) Epidemiology (2067) Evolutionary Biology (24367) Genetics (15621) Genomics (22545) Immunology (17764) Microbiology (40475) Molecular Biology (17208) Neuroscience (88744) Paleontology (667) Pathology (2842) Pharmacology and Toxicology (4834) Physiology (7659) Plant Biology (15175) Scientific Communication and Education (2047) Synthetic Biology (4304) Systems Biology (9834) Zoology (2272)