Front. Comput. Neurosci.Frontiers in Computational NeuroscienceFront. Comput. Neurosci.1662-5188Frontiers Media S.A.10.3389/fncom.2014.00164NeuroscienceOriginal Research ArticleModel reduction of strong-weak neuronsDuBosenSorensenDannyCoxSteven J.*Department of Computational and Applied Mathematics, Rice UniversityHouston, TX, USA
Edited by: Sergey M. Korogod, National Academy of Sciences of Ukraine, Ukraine
Reviewed by: Bartlett W. Mel, University of Southern California, USA; Kresimir Josic, University of Houston, USA
*Correspondence: Steven J. Cox, Department of Computational and Applied Mathematics, Rice University, 6100 Main Street MS 134, Houston, TX 77005, USA e-mail: cox@rice.edu
This article was submitted to the journal Frontiers in Computational Neuroscience.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
We consider neurons with large dendritic trees that are weakly excitable in the sense that back propagating action potentials are severly attenuated as they travel from the small, strongly excitable, spike initiation zone. In previous work we have shown that the computational size of weakly excitable cell models may be reduced by two or more orders of magnitude, and that the size of strongly excitable models may be reduced by at least one order of magnitude, without sacrificing the spatio–temporal nature of its inputs (in the sense we reproduce the cell's precise mapping of inputs to outputs). We combine the best of these two strategies via a predictor-corrector decomposition scheme and achieve a drastically reduced highly accurate model of a caricature of the neuron responsible for collision detection in the locust.
Since Hodgkin and Huxley the neuroscience community has built mathematical models of cells, junctions and circuits as means to both synthesize existing knowledge and to drive further experiments. The complexity of both individual neurons and the networks in which they function has posed serious challenges to those in search of minimal models. The goal of neuronal model reduction is to arrive at a compact description of the cell's “function” and an efficient means of computing its response to physiological stimuli. This is typically accomplished by discovering a smaller equivalent dynamical system and discerning from this a smaller equivalent electrical circuit. See Brunel et al. (2014), Jadi et al. (2014) and Hedrick and Cox (2014) for recent surveys.
We continue our focus, on reduced single cell models that preserve the spatio-temporal structure of their inputs, by providing a detailed synthesis of the active reduction strategy of Kellems et al. (2010) with the quasi-active reduction strategy of Hedrick and Cox (2013). The synthesis is achieved via an elegant method of Rempe and Chopp (2006) for decoupling portions of complex cells and is applied to a caricature of the Lobula Giant Movement Detector (LGMD), the neuron, Peron et al. (2009), responsible for collision detection in the locust. The LGMD has a large, non-spiking dendritic tree that integrates visual input in a retinotopic fashion and funnels this signal to a well defined Spike Initiation Zone (SIZ). Although the structural morphology of the LGMD, and its inputs, has been carefully mapped it is not yet understood what distribution of active and passive conductances permits the cell to discern threatening from, seemingly similar, innocuous visual stimuli. It is hoped that a reduced model will constrain the large parameter space and accelerate the search through this space, and that it will lead to a compact description of the complex task of collision detection as implemented by the full LGMD. For a thorough investigation of the notion of weak excitability in the context of hippocampal pryamidal cells see Golding et al. (2001).
We build and test a detailed (879 compartments) model of the LGMD in §2.1, decouple its branches in §2.2, reduce its active branch in §2.3 and then its quasi-active branches in §2.4. We recouple these two small (3 dimensional) systems in §2.5 and in §3 demonstrate that the drastically reduced system retains the full integrative qualities of the original 879-dimensional model while running 20 times faster.
2. Materials and methods
The caricature of the LGMD neuron raised by Peron et al. (2009) is the rake depicted in Figure 1A. We have numbered its 22 branches and marked its SIZ, in black, near the center of the handle (branch 21) and the joint, in red, where the deck (branch 22) meets the handle. We have chosen a compartment (spatial step) size of dxj = 10 μm and so arrive at a base system with 879 compartments. These are illustrated in Figure 1A and their spatial dimensions are best seen in Figure 1B. We distribute standard sodium, potassium and chloride channels throughout the rake in such a fashion that the tines, branches 1 through 20, weakly integrate synaptic input, funnel it to the deck which then delivers it via the joint to a strongly excitable handle.
A strong-weak neuron. (A) The layout of the rake with its 22 numbered branches, divided into compartments of length dxj = 10 μm. (B) The rest potential, v, obtained by solving Equation 2 subject to the parameters specified in Equations 3 and 4. (C) The response of the full model, Equation 1, at the SIZ (left end of black region in A) and joint (red square in A) to the coherent synaptic input of Equation 5. (D) The response of the full model, Equation 1, at the SIZ (left end of black region in A) and joint (red square in A) to the random synaptic input of Equation 6.
After specifying the full model we decompose it via a predictor-corrector scheme and then apply distinct reduction strategies to the strong and weak parts. Throughout we have used a time step of dt = 0.005 ms.
2.1. The full model
With regard to the rake depicted in Figure 1A, we suppose that the radius of the jth branch is aj = aj(x), where x denotes distance along the branch, and that its associated transmembrane potential is vj = vj(x, t). If the branch contains sodium, potassium and chloride ion channels and is subject to direct current stimulation then Kirchhoff's current law reads
where similar gating equations hold for hj and nj. In addition, we solve Equation 1 subject to sealed ends, current balance at branch points and initial conditions vj(x, 0) = vj(x) and wj(x, 0) = w∞(vj(x)) where vj(x) is the associated rest potential, obtained by solving
will render the tines, branches 1–20, and the deck, branch 22, weakly excitable, while setting
gNa,21(x)={216mS/cm2,200≤x<260μm12mS/cm2,otherwise. and gCl,21=0.3mS/cm2
will make the handle, branch 21, strongly excitable. We have illustrated the resulting rest potential, v, in Figure 1B. We see that the non-uniformity in Equation 4 leads to a depolarized handle and a non-uniform rest potential throughout the remainder of the rake.
We will solve this full system, Equation 1, for two classes of inputs. For the first class, deemed coherent, we simulataneously inject 4 nano-Amperes of current at the midpoint of each tine for nine tenths of a millisecond. In symbols
Istim,j(x,t)=0.004δ(x−200)χ[0.1,1](t),1≤j≤20,
where χ[a, b](t) equals one if a ≤ t ≤ b and equals zero otherwise. For the second class, deemed random, we inject 4 nano-Amperes of current at a random location, and at a random time, on each tine for nine tenths of a millisecond. In symbols
Istim,j(x,t)=0.004δ(x−xj)χ[tj,tj+0.9](t),1≤j≤20,
where the mean of xj is 200 μm and the mean of tj is 5 ms.
In response to coherent stimulus, Equation 5, we see in Figure 1C steady and significant (20 mV) depolarization at the joint (red trace) that is sufficient to drive the handle to spike (blue trace at SIZ). This spike travels down the handle and leads to the second, smaller, depolarization at the joint. The random stimulus, Equation 6, delivers the same amount of current to the rake but spread over space and time. The response at the joint, red trace in Figure 1D, indicate ≈3 mV depolarizations to individual current steps. These are not coherent enough to accumulate in a fashion sufficient to drive the handle to spike. Instead the response at the SIZ, blue trace in Figure 1D, is a filtered attenuated version of the joint trace.
2.2. Branch decomposition
Rempe and Chopp (2006) introduced a rational scheme for decomposing large cells into smaller (typically single branch) regions. They were motivated by the fact that as an action potential travels through a cell, branches on either side of the action potential are relatively quiet and so need not be simulated/computed. As such they devised branch-wise activity measures, in both Rempe and Chopp (2006) and Rempe et al. (2008), that allowed them to build a spatially adaptive numerical scheme that focused resources solely on active branches. One significant advantage of their decomposition is that it permits simultaneous/parallel updating of the active branches. This feature has been successfully exploited by Kozloski and Wagner (2011). Our use of Rempe and Chopp (2006) is however, quite different. For we use their scheme to partition the cell into strong and weak zones that may then be reduced by strategies specific to the dynamics consistent with such zones.
Rempe and Chopp (2006) decompose the cell by giving special attention to those compartments, deemed nodes, at which branches meet. We have illustrated this decomposition on our rake in Figure 2. This spatial decomposition is only useful when coupled with a scheme for properly updating the components in time. Rempe and Chopp (2006) sketch a method that
uses the present branch and node potentials to predict the future node potentials,
updates the branch potentials based on the predicted node potentials,
corrects the node potentials based on the updated branch potentials.
Branch compartment and node labeling to facilitate decoupling via a predictor-corrector scheme. The nodes are colored blue and their potentials are w1 through w21. They occur at the ends of 21 respective branches. The potential in compartment k of branch j is denoted vj(k).
As the success of our method hinges on this predictor-corrector scheme we present it here in some detail.
We distinguish between branches 1 through 21, which are adjacent to a single node, and branch 22, which is adjacent to many. Given the branch, vj(k, t), and node, wj(t), potentials and gating variables at time t we advance the gating variables via the explicit in v implicit in m step
Following Rempe and Chopp (2006) we decouple the node and branch equations in time by making a crude prediction of the nodal potentials by replacing the backward Euler step Equation 8 with the forward Euler step
where wj*(t) denotes our crude prediction of wj(t + dt). We note that Equation 9 may be solved explicitly for
where Γ*j(t + dt) = Γj(40, t + dt) and γ*j(t + dt) = γj(40, t + dt). We then use these predicted nodal potentials to drive the branch updates via
where cj is the node branch coupling vector and Bj is the branch tridiagonal matrix. For the branches adjacent to a single node, j < 22, we find that cj is zero at each compartment except for
This differs from the previous Bj in the sense that it has no free ends (hence 2 terms on the end diagonals) and meets the 20 tines (and hence three terms on those diagonals). The associated coupling term is then zero except at
c22(pj,t)=Gj+∗wj∗(t),j=1,…,21.
Upon updating all branches we may then return to correcting the nodal potentials, now via
With this we may now offer a precise specification of the Predictor-Corrector Algorithm
Given the branch potentials, node potentials and gating variables at time t update the gating variables per Equation 7.
Predict the new values of the node potentials via Equation 10.
Update the branch potentials via Equation 11.
Correct the node potentials via Equation 15. Return to step [1].
2.3. Reduction of the strong part
Following Kellems et al. (2010) we reduce the dynamics in the strong zone, (v21, m21, h21, n21) by the method of Proper Orthogonal Decomposition by collecting snapshots of the membrane potential and associated active current
V=[v21(0)v21(dt)⋯v21(Tfin)] and F=[Iact(0)Iact(dt)⋯Iact(Tfin)]
under a stimulus regime that generates a spike on branch 21. The major features of spike generation and propagation are purportedly captured in the first few singular vectors of V and F. Accordingly we compute the respective singular value decompositions
V=UΣAT and F=WΛCT,
where the matrices of singular vectors, U, A, W and C, are orthonormal and the matrices of singular values, Σ and Λ, are diagonal – and ordered in a decreasing manner.
Our first stab at reduction is to suppose that v21 is well approximated by the first κ columns of U, i.e.,
v21(t)≈Uκv^21(t)
where Uκ denotes the first κ columns of U (from Equation 16) and so the reduced state v^21(t) ∈ ℝκ. On placing this guess in the (spatially discretized version) of Equation 1 we find that the reduced state, v^21, must obey
This provides a clean reduction of the linear spatial coupling between compartments, in the sense that
B˜21≡UκTB21Uκ
is merely κ-by-κ. The non-linearities however are still computed on the full dimensional vector Uκv^21. To address this we distil from Wκ, the first κ columns of W (from Equation 16), κ places along the handle at which it suffices to evaluate the non-linear gating functionals. These places are selected by Discrete Empirical Interpolation as those places at which the singular vectors of F have the greatest content. In particular,
z1=argmax|Wκ(:,1)|P=ez1for i=2:κs=(PTWκ(:,1:i−1))\PTWκ(:,i)r=Wκ(:,i)−Wκ(:,1:i−1)szi=argmax|r|P=[Pezi] end
where ek denotes the kth column of the identity matrix on ℝ39. With these κ places, z = [z1, …, zκ] and their associated permutation matrix P we reduce the gating variables via
m21(t)≈Pm^21(t).h21(t)≈Ph^21(t) and n21(t)≈Pn^21(t)
where gNa(z) denotes the evaluation of gNa, 21 at the compartments indexed by z. As both
R=UκTWκ(PTWκ)−1 and Z=PTUκ
are κ-by-κ we have arrived at a κ-dimesional reduction of the original 39-dimensional active handle. We solve Equation 18, subject to the initial conditions v^21(:, 0) = UTκv and m^21(:, 0) = m∞(Zv^21(:, 0)), via the standard explicit-implicit Euler method
The c21 term in Equation 19 remains the contribution from the nodal potential, w21. Before specifying this we discuss how to reduce the remainder of the branches.
2.4. Reduction of the weak part
In order to perform a single reduction on the remaining branches it is most convenient to gather the variables in the 800 tine compartments and 39 deck compartments into four long vectors v, m, h and n. The first step is then to linearize the full system, Equation 1, about its rest state. More precisely, assuming the stimulus to be order ε we develop the voltage and gating variables
v(x,t)=v¯(x)+εv˜(x,t)+O(ε2) and m(x,t)=m∞(v¯(x))+εm˜(x,t)+O(ε2).
On substituting Equation 20 into Equation 1 and identifying terms of order ε, we find that the so-called quasi-active variables, v˜, m˜, h˜ and n˜ must solve
subject to current balance where the tines meet the deck and to the initial conditions v˜(x, 0) = m˜(x, 0) = 0. On stacking the quasi-active variables in
y=[m˜;h˜;n˜;v˜]
and the stimuli and coupling vector in
u=[Istim,1;Istim,2;⋯;Istim,20;c22]
we may write Equation 21 as an (839 · 4)–dimensional system for y,
y′(t)=Qy(t)+Bu(t),y(0)=0
where the non-zero blocks of Q and B are
Q=(Dm,1Dm,2Dh,1Dh,2Dn,1Dn,2Dv,1Dv,2Dv,3H) and B=(I/Cm).
See §9.4 of Gabbiani and Cox (2010) for the diagonal D matrices and the Hines matrix, H. Given the geometry of the rake we are really only interested in the potential at the joint (recall the red square in Figure 1A). As this is compartment number 820 in the natural ordering, out of the large system Equation 23 we ask only for a good approximation to the joint potential
J(t)=e839·3+820Ty,
where en is the unit vector with a one in element n. Following Hedrick and Cox (2013) we suppose y ≈ ŷ and choose an , with only 4κ columns, that returns an accurate approximation of J. This is done by matching the first κ moments of the full and reduced transfer functions, or, equivalently, via the Arnoldi scheme
x=H\e820X=x/‖x‖for i=1:κ−1x=H\X(:,i) for j=1:ix=x−(X(:,j)Tx)X(:,j) endX=[Xx/‖x‖]end
where, for simplicity, we have chosen our reduced dimension, κ, to agree with that used to reduce the cell's strong zone. We then arrive at the full reducer by tiling this X, i.e.,
On inserting y = ŷ into Equation 23 and using T = I we find that ŷ must obeys
This we solve by backward Euler
and then read off the approximate joint potential via J^(t) = eT839·3+820ŷ(t).
2.5. The reduced strong-weak neuron
It remains only to specify the predictor and corrector updates of the single nodal potential, w21, and to clarify their roles in the coupling vector c21 appearing in the strong reduction, Equation 19, and the coupling vector c22 in the weak reduction, Equation 26, via its presence in the u of Equation 22.
Regarding the expressions for w*21 and w21 in Equation 10 and 15 we note that the required adjacent potentials are readily derived from our independent reductions,
v21(39,t)≈Uκ(39,:)v^(t) and v22(20,t)≈v¯22(20)+J^(t).
The coupling vector c21 is all zero except c21(39, t) = G*21−w*21(t), while the coupling vector c22 is all zero except c22(20, t) = G*21+(w*21(t) − v22(20)).
3. Results
We present in Figure 3 structural components of the strong (A) and weak (B) reductions and in panels (C) and (D) contrast the responses of the full and reduced models, at SIZ and Joint, to the respective coherent and random stimuli used in Figure 1.
The strong-weak reduction of the rake. (A) The singular values of the voltage (V) and active current (F) snapshots. We see that both have decreased by two orders of magnitude by their third index. (B) An illustration of the three columns of the reducer, X, computed in Equation 24, indicating how the reduced system processes the true inputs. (C) Contrasting the response of the full (solid) and reduced (κ = 3) models at the SIZ (black) and at the Joint (red) to identical coherent stimuli. (D) Contrasting the response of the full (solid) and reduced (κ = 3) models at the SIZ (black) and at the Joint (red) to identical random stimuli.
These results were robust to changes in the stimuli that generated the snapshots and to changes in the random stimuli, Equation 6. The reduced system consistently ran in less than 1/20th of the time required by the original system. With κ = 3 the Discrete Empirical Interpolation method identified z = [1, 21, 26] as the compartments along the handle at which to evaluate the gating variables. On comparison to the sodium channel distribution in Equation 4 we note that compartments 21 and 26 are the extent of the SIZ. Regarding Figure 3B we interpret the columns of X as the dendritic filter, seen at the joint, of the true inputs. For columns of X and biophysical interpretations of the elements of TQ see Hedrick and Cox (2013). The errors reported in Figures 3C,D are quite small relative to the original signals in Figures 1C,D and, regarding the timing of critical SIZ events, produce negligible (< 0.1 ms) errors. This then permits us to replace the complex 879 compartment model of Figure 1A with the 8 compartment model of Figure 4.
A shematic of the reduced rake. The true inputs, ui, into the 800 tine compartments are weighted by the columns of X and then summed as they enter the 3 linear nodes, Li, of the reduced weak zone. The linear nodes are fully coupled and a linear combination of their responses contributes to the joint potential, vJ. The 3 fully coupled non-linear nodes, Ni, of the reduced strong zone contribute to both vJ and the SIZ potential vS.
4. Discussion
We have developed and demonstrated a strategy for the systematic reduction of models of strong-weak neurons. In particular, we have replaced a sensory neuron of dimension 879 with a 3-dimensional strong system coupled, via a single node, to a 3-dimensional weak system, Figure 4, and found negligible absolute differences in their voltage responses to complex spatio-temporal inputs while running 20 times faster than the original. We have achieved the strong-weak distinction through significant non-uniformity in the density of sodium channels. This was merely a matter of convenience. The effect can be achieved, see Golding et al. (2001) and Migliore and Shepherd (2002), by a large class of non-uniformities.
The critical assumption permitting our significant reduction is that the bulk of the neuron is weakly excitable - and this means that its response is well approximated by a quasi-active model. The delineation of such systems is of course wrapped up in the equally vexing questions of spike initiation and propagation. For neurons whose dendrites are not sufficiently weak to meet our definition we may apply our strong reduction to each branch and then invoke the activity measures of Rempe and Chopp (2006) and Rempe et al. (2008) to update these branches only when necessary. Regarding scope, our methods are equally suited to synaptic inputs modeled as conductance changes onto trees with arbitrary branching patterns and arbitrary non-uniform channel distributions, see Kellems et al. (2010) and Hedrick and Cox (2013), as well as to inputs via gap junctions, see Hedrick and Cox (2014). These methods can also be readily adapted to incorporate non-uniform distributions of spines and NMDA receptors as well as interaction with the cell's calcium handling machinery of channels, buffers, receptors and pumps.
Conflict of interest statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
This work was partially supported by NSF Grants DMS-1122455 and CCF-1320866 and AFOSR Grant FA9550-12-1-0155.
ReferencesBrunelN.HakimV.RichardsonM. J. E. (2014). Single neuron dynamics and computation. Curr. Opin. Neurobiol. 25, 149–155. 10.1016/j.conb.2014.01.00524492069GabbianiF.CoxS. J. (2010). Mathematics for Neuroscientists. Amsterdam: Academic Press.GoldingN. L.KathW. L.SprustonN. J. (2001). Dichotomy of action-potential backpropagation in CA1 pyramidal neuron dendrites. Neurophysiology86, 2998-3010. Available on line at: http://jn.physiology.org/content/86/6/299811731556HedrickK. R.CoxS. J. (2013). Structure-preserving model reduction of passive and quasi-active neurons. J. Comput. Neurosci. 34, 1–26. 10.1007/s10827-012-0403-y22714391HedrickK. R.CoxS. J. (2014). Morphological reduction of dendritic neurons, in The Computing Dendrite, eds CuntzH.RemmeM.Torben-NielsenB. (New York, NY: Springer), 483–506.JadiM. P.BehabadiB. F.Poleg-PolskyA.SchillerJ.MelB. W. (2014). An augmented two-layer model captures nonlinear analog spatial integration effects in pyramidal neuron dendrites. Proc. IEEE102. 782–798. 10.1109/JPROC.2014.2312671KellemsA. R.ChaturantabutS.SorensenD. C.CoxS. J. (2010). Morphologically accurate reduced order modeling of spiking neurons. J. Comput. Neurosci. 28, 477–494. 10.1007/s10827-010-0229-420300957KozloskiJ.WagnerJ. (2011). An ultrascalable solution to large-scale neural tissue simulation. Front. Neuroinform. 5, 1–21. 10.3389/fninf.2011.0001521954383MiglioreM.ShepherdG. M. (2002). Emerging rules for the distributions of active dendritic conductances. Nat. Rev. Neurosci. 3, 362–370. 10.1038/nrn81011988775PeronS. P.JonesP. W.GabbianiF. (2009). Precise subcellular input retinotopy and its computational consequences in an identified visual interneuron. Neuron63, 830–842. 10.1016/j.neuron.2009.09.01019778511RempeM. J.ChoppD. L. (2006). A predictor–corrector algorithm for reaction-diffusion equations associated with neural activity on branched structures. SIAM J. Sci. Comput. 28, 2139–2161. 10.1137/050643210RempeM. J.SprustonN.KathW. L.ChoppD. L. (2008). Compartmental neural simulations with spatial adaptivity. J. Comput. Neurosci. 25, 465–480. 10.1007/s10827-008-0089-318459041