Model reduction of strong-weak neurons

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.


INTRODUCTION
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.

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 dx j = 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.
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.

THE FULL MODEL
With regard to the rake depicted in Figure 1A, we suppose that the radius of the jth branch is a j = a j (x), where x denotes distance along the branch, and that its associated transmembrane potential is v j = v j (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 h j and n j . In addition, we solve Equation 1 subject to sealed ends, current balance at branch points and initial conditions v j (x, 0) = v j (x) and is the associated rest potential, obtained by solving again subject to sealed ends and current balance at branch points. We concentrate throughout on a single set of parameters. The choice C m = 1.5 μF/cm 2 , R a = 0.05 k cm, a j = 5 μm E Na = 56, E K = −77, E Cl = −68 mV g Na,j = 2 g K,j = 3.6, g Cl,j = 0.9 mS/cm 2 (3) will render the tines, branches 1-20, and the deck, branch 22, weakly excitable, while setting g Na,21 (x) = 216 mS/cm 2 , 200 ≤ x < 260 μm 12 mS/cm 2 , otherwise. and g Cl,21 = 0.3 mS/cm 2 (4) 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 (5) 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 where the mean of x j is 200 μm and the mean of t j 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. 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 (1) uses the present branch and node potentials to predict the future node potentials, (2) updates the branch potentials based on the predicted node potentials, (3) corrects the node potentials based on the updated branch potentials.

BRANCH DECOMPOSITION
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, v j (k, t), and node, w j (t), potentials and gating variables at time t we advance the gating variables via the explicit in v implicit in m step FIGURE 2 | Branch compartment and node labeling to facilitate decoupling via a predictor-corrector scheme. The nodes are colored blue and their potentials are w 1 through w 21 . They occur at the ends of 21 respective branches. The potential in compartment k of branch j is denoted v j (k).
and collect these into We next use these to take a backward Euler step of the associated voltage equation, Equation 1, While at the ends, G j (1 − ) = 0 and With these we may now make sense of the node equation 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 w * j (t) denotes our crude prediction of w j (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 c j is the node branch coupling vector and B j is the branch tridiagonal matrix. For the branches adjacent to a single node, j < 22, we find that c j is zero at each compartment except for while B j is the tridiagonal matrix Turning to the deck, B 22 is the tridiagonal matrix This differs from the previous B j 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 Upon updating all branches we may then return to correcting the nodal potentials, now via which we solve explicitly for 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 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 v 21 is well approximated by the first κ columns of U, i.e., where U κ denotes the first κ columns of U (from Equation 16) and so the reduced statev 21 (t) ∈ R κ . 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 is merely κ-by-κ. The non-linearities however are still computed on the full dimensional vector U κv21 . 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, where g Na (z) denotes the evaluation of g Na,21 at the compartments indexed by z. As both are κ-by-κ we have arrived at a κ-dimesional reduction of the original 39-dimensional active handle. We solve Equation 18, subject to the initial conditionsv 21 ( :, 0) = U T κ v andm 21 ( :, 0) = m ∞ (Zv 21 ( :, 0)), via the standard explicit-implicit Euler method where 21 = R diag(g Na (z).m 3 21 ( :, t + dt).ĥ 21 ( :, t + dt) + g K (z).n 4 21 ( :, t + dt))Z + U T κ diag(g Cl,21 )U κ γ 21 = R (g Na (z).m 3 21 ( :, t + dt).ĥ 21 ( :, t + dt)E Na + g K (z).n 4 21 ( :, t + dt)E K ) + E Cl U T κ g Cl,21 The c 21 term in Equation 19 remains the contribution from the nodal potential, w 21 . Before specifying this we discuss how to reduce the remainder of the branches.

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 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ñ must solve we may write Equation 21 as an (839 · 4)-dimensional system for y, where the non-zero blocks of Q and B are 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 where e n is the unit vector with a one in element n. Following Hedrick and Cox (2013) we suppose y ≈ Xŷ and choose an X , 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\e 820 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 = Xŷ into Equation 23 and using X T X = I we find thatŷ must obeyŝ This we solve by backward Euler and then read off the approximate joint potential viaĴ(t) = e T 839·3+820 Xŷ(t).

THE REDUCED STRONG-WEAK NEURON
It remains only to specify the predictor and corrector updates of the single nodal potential, w 21 , and to clarify their roles in the coupling vector c 21 appearing in the strong reduction, Equation  (20)).

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. 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/20 th 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 X T QX see Hedrick and Cox (2013)  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

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 3dimensional 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 nonuniformity 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.