Mathematical study of a nonlinear neuron model with active dendrites

Abstract: In this work, we have studied an extended version of the cable equation that includes both active and passive membrane properties, under the so-called sealed-end boundary condition. We have thus proved the existence and uniqueness of the weak solution, and defined a novel mathematical form of the somatic cable equation. In particular, we have manipulated the equation set to demonstrate that the diffusion term in the somatic equation is equivalent to the first-order space derivative of the membrane potential in the proximal dendrites. Our conclusion therefore clues how the somatic potential depends on the dynamic of the proximal dendritic segments, and provides the basis for the morphological reduction of neurons without any significant loss of computational properties.


Introduction
The neurons are the computational units of the brain, wherein they form very intricate networks by making over trillion of synaptic connections [1]. Notoriously, the internal state of neuron is encoded as membrane potential, whereas the action potential (i.e. spike) is an atomic event characterized by fast rise and decay of membrane potential, both occurring in 5-10 ms. In 1952, Hogdkin and Huxley defined the first mathematical model of the electric properties of neuron membrane [2] that has brought math and computational modeling closer to neurobiology. By pursuing this integrative approach, it was unraveled the highly non-linear and multi-compartmental nature of the single-neuron computation [3,4], which integrates synaptically-mediated inputs that results in specific patterns of somatic activity [1,5]. In particular, the diffusion of both supra-and sub-threshold signals throughout the entire neuron are substantially modulated by geometric and passive membrane properties of neurite [5][6][7][8], as described by the extended versions of the Hogdkin-Huxley model, which is also known as cable equation [5]. Therefore, these basic neuronal mechanisms altogether implement the logical constructs [1,3,4,8] for brain computation that overall results in consciousness [9,10].
The importance of the electric propagation over the neuron membrane is crucial for inter-neuronal interactions. For synaptic interaction, the action potential propagation to the axon triggers the neurotransmitter release so changing the membrane potential of the post-synaptic neuron. This in turn leads to excitatory or inhibitory interactions that underlay several basic operations of the brain circuits such as synchronization, amplification, and differentiation of the neuronal representation [1].
Also, the electric signals traveling over the neurite diffuse throughout the extra-cellular space, generating an electromagnetic flux that gives rise to emphatic interactions between neurons [1,[11][12][13]. Taking into account the rapid decay of the extra-cellular potential with the distance from neuritic source (< 150 µm) [12], we may speculate that the diffusion over the neuron membrane remains the preponderant factor in realizing the emphatic interactions. Although the importance of these interactions has been underestimated, it was proved that they can readily synchronize entire neuron sets [11][12][13][14], and thus substantially contribute to generate the brain rhythms that are typically recorded during conscious brain processing [15,16]. Therefore, these examples assess the importance of the electric conduction of the neuronal signal over the membrane in the context of large brain circuits.
Moreover, we note that the fact that the axonal activity necessarily follows the somatic response assesses the centrality of the perisomatic region in driving the single-neuron computation. As such, the mathematical characterization of its dynamic is necessary for achieving a better conceptualization of the overall neuronal dynamic and its functional aspects. Indeed, these mechanisms may be artificially mimicked for novel high-tech applications such as memristor, neuromorphic chips, and neural networks.
In this work, we have studied an extended version of the cable equation [2] that includes both active and passive membrane properties, under the so-called sealed-end boundary condition. We have proved the existence and uniqueness of the weak solution, and defined a novel mathematical form of the somatic cable equation. We have then manipulated the equation set to demonstrate that the diffusion term in the somatic equation is equivalent to the first-order space derivative of the membrane potential in the proximal dendrites. Our conclusion therefore clues how the somatic potential depends on the dynamic of the proximal dendritic segments, and provides the basis for the morphological reduction of neurons without significant loss of computational properties.

The biophysical model of the membrane properties
In this section, we have summarized the salient properties of the neuron membrane as described by the cable equation and Hogdkin-Huxley model [2,5], which have been subject to our mathematical manipulation as described in the next sections. Therefore, here we do not provide a comprehensive description which can be accessed through the typical neuroscience literature [5].

The cable equation
In biophysical terms, variation of membrane potential correlates with changes in the selective permeability of the neuron membrane to different ion species, which induces a trans-membrane current flow through the ion channels. In particular, the Hogdkin-Huxley model describes the co-dependence between electrical permittivity of the ion channels and membrane potential in a small portion of neuron membrane [2] that is conceptualized through a characteristic RC circuit ( Figure 1B,  1C). Together, the Hogdkin-Huxley model was extended for explaining the current flow and signal propagation over multiple neurite, as described by the cable equation. Geometrically, each neurite is approximated as a cylinder (with radius a and length L) encased in the neuron membrane (Fig. 1B), which is a bilayered phospholipid barrier endowed by different proteins called ion channels. The ion channel allows the selective passage of ion species through the neuron membrane, inducing a current flow along the longitudinal axis of the neurite, while the phospholipid barrier is electrically insulating. Therefore, the ion channels and phospholipid barrier make the conductive and capacitive properties of the neuron membrane, respectively. Furthermore, the transmembrane current also propagates inside the neurite, along the axial direction, where the internal volume is filled by the cytoplasmic reticulum that has resistive properties. All these properties are modeled by the cable equation.
The cable equation is a partial differential equation defined as where x is the spatial coordinate along the axial direction of the neurite, and indicates the position of a membrane segment.
The term V (in mV) is the membrane potential and is defined as where V i and V e correspond to the values of the electric potential inside and outside of neuron, respectively, while V rest is the resting value of membrane potential. In particular, the difference V i − V e describes the membrane potential, therefore, the term V indicates the deviation of the membrane potential from its resting value. The term τ (in ms) is the membrane time constant and is defined as where r m (in Ωcm) and c m (in F/cm) quantifies the membrane resistance and capacitance, respectively, while R m (in Ωcm 2 ) and C m (in F/cm 2 ) quantifies the specific membrane resistance and capacitance, respectively. In particular, they equate as and The larger C m , the longer the membrane takes to get fully charged or discharged by current injection. The larger R m , the more difficult it is for the current to change V. Therefore, the values of τ provides a measure of the velocity of the electric propagation throughout a neurite segment as well as a measure of its intrinsic reactivity. The λ (in cm) represents the electrotonic length constant and is defined as where r i (in Ωcm) is the axial resistance and R i (in Ω/cm) is the specific axial resistance. The value of λ estimate how far the stationary current propagates along the axial direction of the neurite. Indeed, the higher the R i , the more difficult is for current to propagate along the cytoplasmic medium. The terms I Na and I K are the sodium and potassium currents of the respective ion species flowing across the neuron membrane (Fig. 1C).
The equation 2.1 can be extended as 2.4 The I Na and I K are defined as and where E Na and K are the specific reversal potentials as defined by the Nernst law, g Na and g K are the conductances associated to the channels, m, h, and n describe the states of the gating variables of the channels. The gating particles are distinct in activation or inactivation particles, which increase and decrease the channel permittivity, respectively, during depolarization. In particular, the ion channel underly the action potential generation [5].
To solve the cable equation, one needs to define the boundary condition of the spatial domain of the variable x, which ranges in the interval [0, L]. In literature, they were suggested different types of boundary conditions but only two are the most widely accepted in physiological terms. The most accepted boundary condition is sealed-end boundary conditions, in which one assumes that In order to analyze the voltage clamp experimental setting, where the membrane voltage is maintained at a set level, the boundary condition becomes, for example at x = 0, V(0, t) = V 0 . Now the value V 0 is the membrane potential desired value (in the experiment it is necessary to add a suitable current to hold V 0 ). Remark. At some endpoints x = L, to involve the activation NMDA receptors, one can assume that there is a nonlinear relation between current and applied membrane potential, where g : R → R is a continuous function fitted experimentally by voltage clamping.

Analysis of the space Hodgkin-Huxley model
In order to consider different ionic species and to analyze a more general case we introduce a general space Hodgkin-Huxley model. Let x ∈ [0, L] the distance along the axis of the cylindrical neurite (axon or dendrite without branching), t the time, u(x, t) the membrane potential at point x and time t, N > 1 the number of the ionic species,ḡ j the conductance of the j−th ionic channel, E j the Nernst potential, m j (x, t), h j (x, t), the gating variables corresponding to the dynamics of j−th ionic channel (we are using the Hodgkin-Huxley theory), J(t) a given function (for example related to an external applied current), s j , q j a non-negative integer, we consider the following system The rate functions α j , β j , a j , b j are nonnegative smooth functions, a general expression is given by a rational functions of exponential in u (see e.g. [2]) where c 1 , c 3 , c 4 , u n are nonnegative constants and C 2 , c 5 are positive constants. In the following we will consider an abstract version of the system (3.1), to this aim, we consider the following function space (see [20] for the notation), where H 1 is the Sobolev space with the usual topology.
; there is T > 0 such that exists one and only one solution of the system (3.1) with rate functions α j , β j , a j , b j as in (3.2) for suitable constants c.
Proof. Following from the works of L. Lamberti [17], and M. Mascagni [18,19] we split the system in the non-linear cable equation and in the system of the gating functions m j , h j , Let m 0 ∈ H such that 0 ≤ m o ≤ 1, and consider the problem where the function v is a given function in the space L 2 (0, T ; H), bounded a.e. for (x, t) ∈ [0, L]×[0, T ], with T > 0, and the functions α, β are H−valued integrable nonnegative bounded functions. Let it is easy to prove that the function is the unique solution of the equation (3.7). Taking into account that α, β ≥ 0, we have .
, β(v) ≤ C K , and after some easy computations, for a suitable constant L K , two different functions v 1 , v 2 and two initial conditions m 0,1 , m 0,2 . We point out that the above analysis can be extended to all gating functions m j , h j . Fixed a function v as above, let m j,v , h j,v the corresponding solutions of the gating variables system (3.6), and set the resulting terms defined in (3.4). We have the following estimates, For two different functions v 1 and v 2 and for two sets of initial conditions (m 1 j,0 , h 1 j,0 ), (m 2 j,0 , h 2 j,0 ) for the gating variables, then for a suitable constant C. Starting from a bounded function v, with b v and d v defined in (3.8), we shall first consider the problem which is a linearization of the nonlinear cable equation of (3.1). From known results on linear diffusion equations (see e.g. [21]), the problem (3.10) has a unique solution u such that with L < 1, for every v and w in S z,R,T . Then, there exists T > 0 such that the operator Q is a contraction in S z,R,T and it is well known that there exists a fixed point (S z,R,T is a complete metric space) u = Q(u).
It is easy to see that u, m j (u), h j (u) are the unique solution of (3.1) for t ∈ [0, T ] satisfying u ∈ S z,R,T , moreover u ∈ W. If u * ∈ W is another solution of (3.1), we must have for a small time t. Let T * > 0 be the largest value for which the previous inequality holds, then u * ∈ S z,R,T * and u = u * . If T * < T it is possible to repeat the same argument starting from the new initial point u * 0 = u(T * ) = u * (T * ), and then get a contradiction.
Remark. The dynamics of the gating variables could be described by the ODE's ensuring the existence, uniqueness and boundedness of the functions m j and h j . We have considered a particular form for F j and G j in (3.6).
Using the properties of the operator Q as in Theorem 1. and adopting an argument similar to the proof of the uniqueness of the solution, it is possible to prove a global existence theorem for every finite time T > 0.

A mathematical model for soma and neurites
Taking into account the Evans' findings [23], we have considered a nonlinear version of the Rall's method for morphological complexity reduction. First, we reduce the entire morphology to the equivalent cylinders [24] that connect to the soma (Figure 2). We note that we can operate the morphological reduction under certain condition of spatial symmetry of the dendritic tree [25,26]. However, we note that this reduction methods assure to not alter the dynamic of the soma. In Figure 2 we summarizes the main transformation in constructing the multi-cylinder model. Indeed, under certain symmetry assumptions (see e.g. [26]), each dendrite tree can be represented by an equivalent cylinder. Let M ≥ 1 be the number of neurites (dendrites or axon) connected to the soma and let L i be the length of the equivalent cylinder of the neurite denoted by the index i = 1, 2, . . . , M. We regard any neurite as a one-dimensional continuum reduced to the interval [0, L i ], where the common point x = 0 yields the soma position and x = L i represents the free end. Taking an arbitrary time T > 0 and letting u i : (0, L i ) × (0, T ) → R be the potentials, the transmission in the neurite i is described by the following nonlinear cable equation, Each cable equation is coupled with a system for the gating variables, (3.16) We fixed the initial condition as u i (x, 0) = u i,0 , x ∈ (0, L i ), i = 1, ..., M. We will consider a simple boundary condition at the soma, at x = 0, which corresponds to the experimental situation of the voltage clamp electrophysiological recordings. With voltage-clamp technique intensive research has been conducted to characterize the kinetic properties of ionic currents. Various versions of the method were used to determine the charge carrier, voltage gating, ligand gating, activation, inactivation, recovery, etc. of individual ionic currents. In our model, since all dendrites connect to the soma, the following boundary condition has to be satisfied which are in turn coincident with the somatic voltage membrane u 0 . In light of the Kirchoff's law, applied at x = 0, we have the conservation equation, u s is the soma potential, where a i are known positive constants related to the conductance os the neurite, while I s (t) is an applied current at the soma. The analysis of the linear case for the multi-cylinder model was developed in [27], and [28] when the u s does not change.
where u is the solution of (3.15) with given m i j , h i j , then w is the solution of the following system Now, from the Theorem 4, exists an unique solution of (3.20) and the set S K,T is uniformly bounded, i.e., there is R > 0, such that for any v ∈ S K,T we have |v i (x, t)| ≤ R K a. in the following way. Given v ∈ W, solve the system of ordinary differential equation (3.16) and find the gating functions m i j (v), h i j (v); finally find the solution u of the nonlinear cable equation (3.15) and put Q(v) = u. From Theorem 3, it is possible to state for a suitable constant K 0 , and from Theorem 4 that there exists K 1 > 0 such that Then Q : S K,T → S K,T for T (K, K1) < (K/K 1 ) 2 . Moreover Q is a contraction for sufficiently small T * (K, K 1 ). Due to the uniform boundedness of S K,T , we can restart from T * until we get a solution for the whole interval [0, T ].
Remark. Note that, by interpolation, we obtain that then the equation for the somatic potential holds for any t ∈ (0, T ). Furthermore, under the appropriate hypotheses for the current I s , one can prove existence and uniqueness of the solution for the equation (3.21) by the classic theory of ordinary differential equations.
Remark. For the axon, the axial current propagates only from the soma, therefore, we can consider again the reaction-diffusion equation defined as where u A is the axon membrane potential, L A is the axon length, under the boundary conditions u(0, t) = u s (t) and ∂u A /∂x(L A , t) = 0. Remark. We analyzed an approximation of the cable equation using classical finite-difference method. The original equation (2.4) defined by Hogdkin and Huxley can be rewritten as (we consider an uniform mesh with space step ∆x), v s (t + ∆t) = v s (t) + ∆t The term v j i , which describes the first compartment of the dendrite, is coincident with v s . Hence the terms v j i is the same value as v j−1 i because the soma can be approximate by an isopotential surface. Therefore, the equation (3.22) can be rewritten as If we assume that which coincides with the finite-difference form of the equation (3.15). Therefore, one implicitly applies the equation (3.15) by implementing the original cable equation (2.4) in form of finite-difference derivative.

A numerical test
We consider a preliminary numerical test where we compare passive/active dendrites tree with different morphology. To carry out our simulations, we used two neuron morphologies, one concerning of the full dendrite tree, another which consists of the reduced morphology by applying the Ralls law (see Figure 3A). Both morphology includes a soma which is approximated as an isopotential sphere of 80 mm in diameter. The full morphology is formed by 12 dendrites approximated as cylinders. There are 4 and 8 dendrites arranged in the first and second order branches, respectively. For both orders, cylinders are 100 mm in length, whereas the diameters are 20 mm 12.6 mm for the first and second order dendrites, respectively. Thus, the diameters matched the 2/3 Ralls power law. The reduced morphology is formed by 4 identical cylinder that are 200 mm in length and 20 mm in diameter. Passive and active membrane properties were identical for both morphology and uniform over the entire neuron, as summarized in Table 1.  We have considered only two species of ionic currents as described in Hodgkin and Huxley's model: sodium current (Na) with total conductance g Na m 3 , and potassium current (K) with conductance g K n 4 . We have also considered two experimental cases for the injected current: a case with a weak current, and a case with a strong current. For numerical tests we used a multi-compartment approach we a standard time integration Ruge-Kutta method. Figure 3B shows the influence of the different dendrite morphology on the dynamics, in particular differences are observed in the hyper-polarization phase of the action potentials, and in the firing rate. This effect is more evident in the case of stronger currents. During a one-second numerical simulation, the morphology with more branches produces 60 spikes in the weak current case, and 80 spikes for the strong one. While the dynamics of the reduced morphology consists of 70 spikes for the weak current and of 100 spikes corresponding to the strong current.

Conclusion
The dynamic underlying the neuronal signal over the membrane has been a subject of great interest [29], since it underlies all the brain functioning based on non-linear and compartmental computation [3,4,8] as well as synaptic and ephaptic interactions [1,[11][12][13][14][15][16]. Starting from the cable equation [5], here we have re-defined it for the soma and the neurites in its proximity, demonstrating that the membrane potential diffuses linearly from the neurites to the soma. In particular, we have demonstrated the existence and uniqueness of the solution for the cable equation as represented by a mathematical system of reaction-diffusion equations under suitable boundary conditions. We have demonstrated the this set of differential equations is well-posed, which is an important property that provide the basis for further extensions. This provides the basis for further extensions of the model.
In particular, in the future, we may consider a network of realistic neuron where the equations for the synaptic terminals are subject to non-linear boundary conditions; we may apply the model to reduce the morphological complexity of neuron without causing loss of the computational properties. This would benefit us over the artificial replication of the biophysical mechanisms but also to analyze the intrinsic properties of single-neurons. For example, in previous works, computational modeling based on the Hogdkin-Huxley equations was applied for analyzing the intrinsic properties of neuron that determine the action potential propagation over neurites [7]. Instead, here we have unraveled a biophysical concept by re-formalizing the cable equation, giving novel insight into how the somatic response depends linearly on the proximal neurites. In the next work, we aim to address when this linearity extends to the non-proximal portions of the dendrites, in order to understand whether there exist specific values of the biophysical properties that may affect the nature of signal propagation rather than just the speed or the dispersion as it is more commonly assumed in literature.