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.

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

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

FIGURE 1

Figure 1. A strong-weak neuron. (A) The layout of the rake with its 22 numbered branches, divided into compartments of length dx_j = 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 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

$\begin{array}{cc}\begin{array}{l}{C}_m{\partial }_t{v}_j=\frac{1}{2a_j{R}_a}{\partial }_x(a_j^2{\partial }_x{v}_j)+I_{stim,j}/(2\pi a_j)\hfill \\ -g_{Na,j}{m}_j^3{h}_j(v_j-E_{Na})-g_{K,j}{n}_j^4(v_j-E_K)\hfill \\ -g_{Cl,j}(v_j-E_{Cl})\hfill \\ {\partial }_t{m}_j=(m_{\infty }(v_j)-m_j)/{\tau }_m(v_j),\hfill \end{array}& (1)\end{array}$

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 w_j(x, 0) = w_∞(v_j(x)) where v_j(x) is the associated rest potential, obtained by solving

$\begin{array}{cc}\begin{array}{l}{(a_j^2(x){{\overline{v}}^{\prime }}_j(x))}^{\prime }=2a_j{R}_a\{g_{Na,j}(x)m_{\infty }^3({\overline{v}}_j(x))h_{\infty }({\overline{v}}_j(x))({\overline{v}}_j(x)-E_{Na})\hfill \\ +g_{K,j}(x)n_{\infty }^4({\overline{v}}_j(x))({\overline{v}}_j(x)-E_K)\hfill \\ +g_{Cl,j}(x)({\overline{v}}_j(x)-E_{Cl})\},\hfill \end{array}& (2)\end{array}$

again subject to sealed ends and current balance at branch points. We concentrate throughout on a single set of parameters. The choice

$\begin{array}{cc}\begin{array}{l}{C}_m=1.5\mu F/c{m}^2,R_a=0.05k\Omega cm,a_j=5\mu m\hfill \\ E_{Na}=56,E_K=-77,E_{Cl}=-68mV\hfill \\ g_{Na,j}=2g_{K,j}=3.6,g_{Cl,j}=0.9mS/c{m}^2\hfill \end{array}& (3)\end{array}$

will render the tines, branches 1–20, and the deck, branch 22, weakly excitable, while setting

$\begin{array}{cc}\begin{array}{l}{g}_{Na,21}(x)=\{\begin{array}{ll}216mS/c{m}^2,\hfill & 200\le x<260\mu m\hfill \\ 12mS/c{m}^2,\hfill & \text{otherwise}.\hfill \end{array}\\ \text{and}{g}_{Cl,21}=0.3mS/c{m}^2\end{array}& (4)\end{array}$

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

$\begin{array}{cc}{I}_{stim,j}(x,t)=0.004\delta (x-200){\chi }_{[0.1,1]}(t),1\le j\le 20,& (5)\end{array}$

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

$\begin{array}{cc}{I}_{stim,j}(x,t)=0.004\delta (x-x_j){\chi }_{[t_j,t_j+0.9]}(t),1\le j\le 20,& (6)\end{array}$

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.

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

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

FIGURE 2

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₁ through w₂₁. They occur at the ends of 21 respective branches. The potential in compartment k of branch j is denoted v_j(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, 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

and collect these into

$\begin{array}{l}{\Gamma }_j(k,t+dt)=g_{Na,j}(k)m_j^3(k,t+dt)h_j(k,t+dt)\hfill \\ +g_{K,j}(k)n_j^4(k,t+dt)+g_{Cl,j}(k)\hfill \\ {\gamma }_j(k,t+dt)=g_{Na,j}(k)m_j^3(k,t+dt)h_j(k,t+dt)E_{Na}\hfill \\ +g_{K,j}(k)n_j^4(k,t+dt)E_K+g_{Cl,j}(k)E_{Cl}.\hfill \end{array}$

We next use these to take a backward Euler step of the associated voltage equation, Equation 1,

$\begin{array}{l}{G}_j(k-)v_j(k-1,t+dt)-(G_j(k-)+G_j(k+))\hfill \\ v_j(k,t+dt)+G_j(k+)v_j(k+1,t+dt)=\hfill \\ \mu (v_j(k,t+dt)-v_j(k,t))+{\Gamma }_j(k,t+dt)v_j(k,t+dt)\hfill \\ -{\gamma }_j(k,t+dt)+I_{stim,j}(k,t+dt)/A_j(k)\hfill \end{array}$

where

$\begin{array}{l}{G}_j(k\pm )=\frac{a_j(k)}{R_{a}d{x}_j^2}\frac{a_j^2(k\pm 1)}{a_j^2(k)+a_j^2(k\pm 1)},\mu =C_m/dt\hfill \\ \text{and}{A}_j(k)=2\pi a_j(k)d{x}_j.\hfill \end{array}$

While at the ends, G_j(1 −) = 0 and

$\begin{array}{l}{G}_j(40-)\equiv G_{j-}^{\ast }=\frac{a_j(40)}{R_{a}d{x}_j^2}\frac{a_j^2(39)}{a_j^2(40)+a_j^2(39)}\\ G_j(40+)\equiv G_{j+}^{\ast }=\frac{a_j(40)}{2R_{a}d{x}_j}\frac{d{x}_{22}^2}{2a_j^2(40)d{x}_j+d{x}_{22}^2{a}_{22}(p_j)}\\ \text{where}{p}_j=\{\begin{array}{ll}2j-1,\hfill & 1\le j\le 20\hfill \\ 20,\hfill & j=21.\hfill \end{array}\end{array}$

With these we may now make sense of the node equation

$\begin{array}{cc}\begin{array}{l}{G}_{j-}^{\ast }{v}_j(39,t+dt)-(G_{j-}^{\ast }+G_{j+}^{\ast })\\ w_j(t+dt)+G_{j+}^{\ast }{v}_{22}(p_j,t+dt)=\\ \mu (w_j(t+dt)-w_j(t))+{\Gamma }_j(40,t+dt)\\ w_j(t+dt)-{\gamma }_j(40,t+dt)\end{array}& (8)\end{array}$

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

$\begin{array}{cc}\begin{array}{l}{G}_{j-}^{\ast }{v}_j(39,t)-(G_{j-}^{\ast }+G_{j+}^{\ast })w_j(t)+G_{j+}^{\ast }{v}_{22}(p_j,t)=\\ \mu (w_j^{\ast }(t)-w_j(t))+{\Gamma }_j(40,t+dt)w_j^{\ast }(t)-{\gamma }_j(40,t+dt)\end{array}& (9)\end{array}$

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

$\begin{array}{cc}{c}_j(39,t)=G_{j-}^{\ast }{w}_j^{\ast }(t),& (12)\end{array}$

while B_j is the tridiagonal matrix

$\begin{array}{l}{B}_j(1,1:2)=[-G_j(1+\text{})G_j(1+\text{})]\hfill \\ B_j(k,k-1:k+1)=[G_j(k-\text{})-(G_j(k-\text{})+G_j(k+\text{}))\hfill \\ G_j(k+\text{})],1<k<39,\hfill \\ B_j(39,38:39)=[G_j(39-\text{})-(G_j(39-\text{})\hfill \\ +G_j(39+))].\hfill \end{array}$

Turning to the deck, B₂₂ is the tridiagonal matrix

$\begin{array}{l}{B}_{22}(1,1:2)=[-(G_{22}(1+\text{})+G_{1+}^{\ast }\text{})G_j(1+\text{})]\hfill \\ {}_{22}(k,k-1:k+1)=[G_{22}(k-\text{})-(G_{22}(k-\text{})\hfill \\ +G_{22}(k+\text{}))G_{22}(k+\text{})],\hfill \\ 1<k<39,k\ne p_j\hfill \\ B_{22}(k,k-1:k+1)=[G_{22}(k-\text{})-(G_{22}(k-\text{})\hfill \\ +G_{22}(k+\text{})+G_{j+}^{\ast })G_{22}(k+\text{})],\hfill \\ 1<k<39,k=p_j\hfill \\ B_{22}(39,38:39)=[G_{22}(39-\text{}\text{})-(G_{22}(39-\text{})+G_{20+}^{\ast })].\hfill \end{array}$

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

$\begin{array}{cc}{c}_{22}(p_j,t)=G_{j+}^{\ast }{w}_j^{\ast }(t),j=1,\dots ,21.& (13)\end{array}$

Upon updating all branches we may then return to correcting the nodal potentials, now via

$\begin{array}{cc}\begin{array}{l}{G}_{j-}^{\ast }{v}_j(39,t+dt)-(G_{j-}^{\ast }+G_{j+}^{\ast })\hfill \\ w_j(t+dt)+G_{j+}^{\ast }{v}_p(p_j,t+dt)=\hfill \\ \mu (w_j(t+dt)-w_j^{\ast }(t))+{\Gamma }_j^{\ast }(t+dt)\hfill \\ w_j(t+dt)-{\gamma }_j^{\ast }(t+dt)\hfill \end{array}\hfill & (14)\end{array}$

which we solve explicitly for

With this we may now offer a precise specification of the Predictor-Corrector Algorithm

[1] Given the branch potentials, node potentials and gating variables at time t update the gating variables per Equation 7.

[2] Predict the new values of the node potentials via Equation 10.

[3] Update the branch potentials via Equation 11.

[4] 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, (v₂₁, m₂₁, h₂₁, n₂₁) by the method of Proper Orthogonal Decomposition by collecting snapshots of the membrane potential and associated active current

$\begin{array}{l}{I}_{act}(t)\equiv g_{Na,21}.m_{21}^3(t).h_{21}(t).(v_{21}(t)-E_{Na})\hfill \\ +g_{K,21}.n_{21}^4(t).(v_{21}(t)-E_K)\hfill \end{array}$

in

$\begin{array}{l}V=[v_{21}(0)v_{21}(dt)\cdots v_{21}(T_{fin})]\text{and}\hfill \\ F=[I_{act}(0)I_{act}(dt)\cdots I_{act}(T_{fin})]\hfill \end{array}$

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

$\begin{array}{cc}V=U\Sigma A^T\text{and}F=W\Lambda C^T,& (16)\end{array}$

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₂₁ is well approximated by the first κ columns of U, i.e.,

$v_{21}(t)\approx U_{\kappa }{\widehat{v}}_{21}(t)$

where U_κ denotes the first κ columns of U (from Equation 16) and so the reduced state $\widehat{v}$₂₁(t) ∈ ℝ^κ. On placing this guess in the (spatially discretized version) of Equation 1 we find that the reduced state, $\widehat{v}$₂₁, must obey

$\begin{array}{cc}\begin{array}{l}{C}_m{{\widehat{v}}^{\prime }}_{21}=U_{\kappa }^T{B}_{21}{U}_{\kappa }{\widehat{v}}_{21}-U_{\kappa }^T\{g_{Na,21}.m_{21}^3.h_{21}.(U_{\kappa }{\widehat{v}}_{21}-E_{Na})\hfill \\ +g_{K,21}.n_{21}^4.(U_{\kappa }{\widehat{v}}_{21}-E_K)+g_{Cl,21}(U_{\kappa }{\widehat{v}}_{21}-E_{Cl})-c_{21}\}\hfill \\ {m^{\prime }}_{21}=(m_{\infty }(U_{\kappa }{\widehat{v}}_{21})-m_{21})./{\tau }_m(U_{\kappa }{\widehat{v}}_{21}).\hfill \end{array}\hfill & (17)\end{array}$

This provides a clean reduction of the linear spatial coupling between compartments, in the sense that

${\tilde{B}}_{21}\equiv U_{\kappa }^T{B}_{21}{U}_{\kappa }$

is merely κ-by-κ. The non-linearities however are still computed on the full dimensional vector U_κ$\widehat{v}$₂₁. 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,

$\begin{array}{l}{z}_1=\text{argmax}|W_{\kappa }(\text{}:,1)|\hfill \\ P=e_{z_1}\hfill \\ \text{for}i=2:\kappa \hfill \\ s=(P^T{W}_{\kappa }(\text{}:,1:i-1))\backslash P^T{W}_{\kappa }(\text{}:,i)\hfill \\ r=W_{\kappa }(\text{}:,i)-W_{\kappa }(\text{}:,1:i-1)s\hfill \\ z_i=\text{argmax}\left|r\right|\hfill \\ P=[P{e}_{z_i}]\hfill \\ \text{end}\hfill \end{array}$

where e_k denotes the kth column of the identity matrix on ℝ³⁹. With these κ places, z = [z₁, …, z_κ] and their associated permutation matrix P we reduce the gating variables via

$m_{21}(t)\approx P{\widehat{m}}_{21}(t).h_{21}(t)\approx P{\widehat{h}}_{21}(t)\text{and}{n}_{21}(t)\approx P{\widehat{n}}_{21}(t)$

and so bring Equation 17 to

$\begin{array}{cc}\begin{array}{l}{C}_m{{\widehat{v}}^{\prime }}_{21}={\tilde{B}}_{21}{\widehat{v}}_{21}-R\{g_{Na}(z).{\widehat{m}}_{21}^3.{\widehat{h}}_{21}.(Z{\widehat{v}}_{21}-E_{Na})\hfill \\ +g_K(z).{\widehat{n}}_{21}^4.(Z{\widehat{v}}_{21}-E_K)\}-U_{\kappa }^T{g}_{Cl,21}(U_{\kappa }{\widehat{v}}_{21}-E_{Cl})\hfill \\ +U_{\kappa }^T{c}_{21}\hfill \\ {{\widehat{m}}^{\prime }}_{21}=(m_{\infty }(Z{\widehat{v}}_{21})-{\widehat{m}}_{21})./{\tau }_m(Z{\widehat{v}}_{21})\hfill \end{array}\hfill & (18)\end{array}$

where g_Na(z) denotes the evaluation of g_{Na, 21} at the compartments indexed by z. As both

$R=U_{\kappa }^T{W}_{\kappa }{(P^T{W}_{\kappa })}^{-1}\text{and}Z=P^T{U}_{\kappa }$

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 $\widehat{v}$₂₁(:, 0) = U^T_κv and $\widehat{m}$₂₁(:, 0) = m_∞(Z$\widehat{v}$₂₁(:, 0)), via the standard explicit-implicit Euler method

where

$\begin{array}{l}{\Gamma }_{21}=R\text{diag}(g_{Na}(z).{\widehat{m}}_{21}^3(\text{}:,t+dt).{\widehat{h}}_{21}(\text{}:,t+dt)\hfill \\ +g_K(z).{\widehat{n}}_{21}^4(\text{}:,t+dt))Z+U_{\kappa }^T\text{diag}(g_{Cl,21})U_{\kappa }\hfill \\ {\gamma }_{21}=R(g_{Na}(z).{\widehat{m}}_{21}^3(\text{}:,t+dt).{\widehat{h}}_{21}(\text{}:,t+dt)E_{Na}\hfill \\ +g_K(z).{\widehat{n}}_{21}^4(\text{}:,t+dt)E_K)+E_{Cl}{U}_{\kappa }^T{g}_{Cl,21}\hfill \end{array}$

The c₂₁ term in Equation 19 remains the contribution from the nodal potential, w₂₁. 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

$\begin{array}{cc}\begin{array}{l}v(x,t)=\overline{v}(x)+\epsilon \tilde{v}(x,t)+O({\epsilon }^2)\hfill \\ \text{and}m(x,t)=m_{\infty }(\overline{v}(x))+\epsilon \tilde{m}(x,t)+O({\epsilon }^2).\hfill \end{array}\hfill & (20)\end{array}$

On substituting Equation 20 into Equation 1 and identifying terms of order ε, we find that the so-called quasi-active variables, $\tilde{v}$, $\tilde{m}$, $\tilde{h}$ and $\tilde{n}$ must solve

$\begin{array}{cc}\begin{array}{l}{C}_m{\partial }_t\tilde{v}=\frac{1}{2R_{a}a}{\partial }_x(a^2{\partial }_x\tilde{v})-(g_{Na,j}{m}_{\infty }^3(\overline{v})h_{\infty }(\overline{v})+g_{K,j}{n}_{\infty }^4(\overline{v})\hfill \\ +g_{Cl,j})\tilde{v}-3g_{Na,j}{m}_{\infty }^2(\overline{v})h_{\infty }(\overline{v})(\overline{v}-E_{Na})\tilde{m}\hfill \\ -g_{Na,j}{m}_{\infty }^3(\overline{v})(\overline{v}-E_{Na})\tilde{h}-4g_{K,j}{n}_{\infty }^3(\overline{v})(\overline{v}-E_K)\tilde{n}\hfill \\ +I_{stim}/(2\pi a)\hfill \\ {\partial }_t\tilde{m}=({m^{\prime }}_{\infty }(\overline{v})-\tilde{m})/{\tau }_m(\overline{v})\hfill \end{array}\hfill & (21)\end{array}$

subject to current balance where the tines meet the deck and to the initial conditions $\tilde{v}$(x, 0) = $\tilde{m}$(x, 0) = 0. On stacking the quasi-active variables in

$y=[\tilde{m};\tilde{h};\tilde{n};\tilde{v}]$

and the stimuli and coupling vector in

$\begin{array}{cc}u=[I_{stim,1};I_{stim,2};\cdots ;I_{stim,20};c_{22}]& (22)\end{array}$

we may write Equation 21 as an (839 · 4)–dimensional system for y,

$\begin{array}{cc}{y}^{\prime }(t)=Qy(t)+Bu(t),y(0)=0& (23)\end{array}$

where the non-zero blocks of Q and B are

$Q=\left(\begin{array}{cccc}{D}_{m,1}& & & D_{m,2}\\ & D_{h,1}& & D_{h,2}\\ & & D_{n,1}& D_{n,2}\\ D_{v,1}& D_{v,2}& D_{v,3}& H\end{array}\right)\text{and}B=\left(\begin{array}{c}\\ \\ \\ I/C_m\end{array}\right).$

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)=e_{839·3+820}^{T}y,$

where e_n 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

$\begin{array}{cc}\begin{array}{l}x=H\backslash e_{820} \hfill \\ X=x/\Vert x\Vert \hfill \\ \text{for}i=1:\kappa -1\hfill \\ x=H\backslash X(\text{}:,i)\hfill \\ \text{for} j=1:i\hfill \\ x=x-(X{(\text{}:,j)}^{T}x)X(\text{}:,j)\hfill \\ \text{end} \hfill \\ X=[Xx/\Vert x\Vert ]\hfill \\ \text{end} \hfill \end{array}\hfill & (24)\end{array}$

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 $\widehat{J}$(t) = e^T_839·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, w₂₁, and to clarify their roles in the coupling vector c₂₁ appearing in the strong reduction, Equation 19, and the coupling vector c₂₂ in the weak reduction, Equation 26, via its presence in the u of Equation 22.

Regarding the expressions for w*₂₁ and w₂₁ in Equation 10 and 15 we note that the required adjacent potentials are readily derived from our independent reductions,

$v_{21}(39,t)\approx U_{\kappa }(39,:)\widehat{v}(t)\text{and}{v}_{22}(20,t)\approx {\overline{v}}_{22}(20)+\widehat{J}(t).$

The coupling vector c₂₁ is all zero except c₂₁(39, t) = G*₂₁₋w*₂₁(t), while the coupling vector c₂₂ is all zero except c₂₂(20, t) = G*₂₁₊(w*₂₁(t) − v₂₂(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.

FIGURE 3

Figure 3. 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/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 ^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.

FIGURE 4

Figure 4. A shematic of the reduced rake. The true inputs, u_i, into the 800 tine compartments are weighted by the columns of X and then summed as they enter the 3 linear nodes, L_i, of the reduced weak zone. The linear nodes are fully coupled and a linear combination of their responses contributes to the joint potential, v_J. The 3 fully coupled non-linear nodes, N_i, of the reduced strong zone contribute to both v_J and the SIZ potential v_S.

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.

Acknowledgments

This work was partially supported by NSF Grants DMS-1122455 and CCF-1320866 and AFOSR Grant FA9550-12-1-0155.

References

Brunel, N., Hakim, V., and Richardson, M. J. E. (2014). Single neuron dynamics and computation. Curr. Opin. Neurobiol. 25, 149–155. doi: 10.1016/j.conb.2014.01.005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gabbiani, F., and Cox, S. J. (2010). Mathematics for Neuroscientists. Amsterdam: Academic Press.

Google Scholar

Golding, N. L., Kath, W. L., and Spruston, N. J. (2001). Dichotomy of action-potential backpropagation in CA1 pyramidal neuron dendrites. Neurophysiology 86, 2998-3010. Available on line at: http://jn.physiology.org/content/86/6/2998

Pubmed Abstract | Pubmed Full Text

Hedrick, K. R., and Cox, S. J. (2013). Structure-preserving model reduction of passive and quasi-active neurons. J. Comput. Neurosci. 34, 1–26. doi: 10.1007/s10827-012-0403-y

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hedrick, K. R., and Cox, S. J. (2014). “Morphological reduction of dendritic neurons,” in The Computing Dendrite, eds H. Cuntz, M. Remme, and B. Torben-Nielsen (New York, NY: Springer), 483–506.

Google Scholar

Jadi, M. P., Behabadi, B. F., Poleg-Polsky, A., Schiller, J., and Mel, B. W. (2014). An augmented two-layer model captures nonlinear analog spatial integration effects in pyramidal neuron dendrites. Proc. IEEE 102. 782–798. doi: 10.1109/JPROC.2014.2312671

CrossRef Full Text | Google Scholar

Kellems, A. R., Chaturantabut, S., Sorensen, D. C., and Cox, S. J. (2010). Morphologically accurate reduced order modeling of spiking neurons. J. Comput. Neurosci. 28, 477–494. doi: 10.1007/s10827-010-0229-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kozloski, J., and Wagner, J. (2011). An ultrascalable solution to large-scale neural tissue simulation. Front. Neuroinform. 5, 1–21. doi: 10.3389/fninf.2011.00015

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Migliore, M., and Shepherd, G. M. (2002). Emerging rules for the distributions of active dendritic conductances. Nat. Rev. Neurosci. 3, 362–370. doi: 10.1038/nrn810

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Peron, S. P., Jones, P. W., and Gabbiani, F. (2009). Precise subcellular input retinotopy and its computational consequences in an identified visual interneuron. Neuron 63, 830–842. doi: 10.1016/j.neuron.2009.09.010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Rempe, M. J., and Chopp, D. L. (2006). A predictor–corrector algorithm for reaction-diffusion equations associated with neural activity on branched structures. SIAM J. Sci. Comput. 28, 2139–2161. doi: 10.1137/050643210

CrossRef Full Text | Google Scholar

Rempe, M. J., Spruston, N., Kath, W. L., and Chopp, D. L. (2008). Compartmental neural simulations with spatial adaptivity. J. Comput. Neurosci. 25, 465–480. doi: 10.1007/s10827-008-0089-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar