An electromechanical model of neuronal dynamics using Hamilton's principle

Introduction

Brain tissue is an inhomogeneous, multi-scale composite material composed of interconnected networks of blood vessels, neuron, and glia cells submerged in cerebrospinal fluid. Effects of mechanical and/or electro-chemical stresses and deformations on brain vary widely depending on the cell types, mechanical and bio-chemical characteristics of the cells, as well as cell's mechanosensitivity and mechanotransduction abilities. For instance, brain damage may take many different forms. For neurons, damage might include breakage of cytoskeleton networks in dendrites or axons, membrane rupture, separation of synaptic connections, or severance of dendritic or axonal projections. For the vascular system, damage might be puncture of macro or micro capillaries, or restrictions that alter perfusion on various scales. Depending on severity, most of these mechanical injuries will be followed by short or long term chemical imbalances and/or functional impairments or even death.

Given the high complexity of brain's structure and dynamics, designing, performing, and interpreting experiments on brain in vivo at various time and length scales continue to be very challenging and as a result the mechanisms that govern the interconnected electro-chemo-mechanical processes that result in brain damage and ultimately loss of functionality remain poorly understood. Mathematical models and corresponding computer simulations can increase our comprehension on brain damage processes (and, in general, on neurological diseases and neurodegeneration) and help design better experiments for measurements and hypothesis testing that ultimately will lead to improved medical diagnostic and therapeutic protocols. In the last few decades a multitude of mathematical models have been proposed to study brain biomechanics and, independently, brain bio-chemistry at cell as well as tissue levels. The majority of these models are mentioned in the comprehensive reviews of brain biomechanics and mechanobiology by Goldsmith (2001) and Goriely et al. (2015). Recently, models that link brain biomechanics to its bio-chemistry have also started to be developed (Drapaca and Fritz, 2012; Lang et al., 2015). Such coupled models are essential in enhancing the understanding of brain mechanisms such as the onset of normal pressure hydrocephalus due to ionic imbalances and in the absence of an elevated intracranial pressure (Drapaca and Fritz, 2012), and the propagation of damage in brain tissue caused by edema and lack of proper oxygenation (Lang et al., 2015). However, the lack of experimental data required by these very advanced mathematical models as well as the complexity of the corresponding computations render these models hard to use in today's clinical applications. In addition, these coupled models have been built at tissue level and thus they cannot predict the mechano-chemical responses of brain cells to mechanical and/or electro-chemical events that happen at tissue and organ scales.

The latest survey of the literature on brain biomechanics and mechanobiology by Goriely et al. (2015) emphasizes the current need in the field of brain research for the development of “bottom-up” mathematical models that link brain mechanics and function at each relevant length scale as well as across scales, incorporate anatomically accurate geometry and connections of cells and cerebral vasculature, and ultimately allow information from molecular and cellular levels to propagate to tissue and organ levels and vice versa. One possible first step in building such a bottom-up model is to start at the cell level and create an electromechanical model of neuronal dynamics. The aim of this paper is therefore to develop a lower-dimensional electromechanical model of a neuron which (1) is simple enough so that its predictions may be experimentally verified, and (2) could be used as a foundation model for more advanced multi-scaling (bottom-up) mathematical models. We assume that the electro-chemical activity of a neuron is described by the classic Hodgkin–Huxley equations (Hodgkin and Huxley, 1952) and that from a mechanical point of view the neuron behaves like a linear visco-elastic Kelvin–Voigt solid. The assumption of linear viscoelastic neuron is supported by experimental evidence reported by Lu et al. (2006) and Grevesse et al. (2015). In order to couple the Kelvin–Voigt mechanical model and the Hodgkin–Huxley electric model we will use a constraint Lagrangian formulation and the non-conservative form of Hamilton's principle. This approach will give us the coupled equations of motion by minimizing a special integral functional (action) whose integrand is made of kinetic and potential energies (Lagrangian) and the work done by the forces acting on the neuron. Although Hamilton's principle has been used in classical mechanics for a very long time (see for instance Lanczos, 1986), and recently has been applied to model neuronal electro-chemical activities (Dickel, 1989; Paninski, 2006; Wilson and Steyn-Ross, 2008; Chuankui, 2012) and ion transport through cell's membrane (Eisenberg et al., 2010), the principle has not been used to link neuron's mechanics and its electro-chemistry until now. The proposed electromechanical model has the following desirable features: (1) incorporates relevant macroscopic (cell level) and microscopic (ionic level) mechanical and electrical information, (2) facilitates the study of the dynamics of neuronal stiffness due to the evolutions of microstructural components, and (3) highlights neuronal mechanotransduction. We test the performance of our model in numerical simulations of neuronal mechanical insults. Although today it is well-known and accepted that traumatic brain injuries change the mechanics and electrophysiology of neurons on short and long time scales (see for instance Goriely et al., 2015, and the references within), the focus of the experimentalists as well as the modelers has been primarily on the mechanical characterization of the neuronal damage, and therefore a direct liaison between the neuronal mechanical properties and its altered functions has not been established yet. Our numerical simulations clearly show neuronal mechanotransduction: for initially applied displacements and speeds of magnitudes comparable to the size of the neuron, action potentials are observed, while very fast initially applied speeds (jabbing) inhibit the action potentials and this case might describe one possible neuronal damage dynamics following a serious mechanical injury. In addition, we notice that our proposed dynamics for the stiffness of a neuron appears to be in agreement with the experimental observations of healthy neurons reported by Zou et al. (2013).

The paper is organized as follows. In Section Mathematical Model we present our mathematical model, and in Section Results we show some relevant numerical simulations. The paper ends with a section containing concluding remarks and future directions.

Mathematical Model

We model the axon as a axi-symmetric circular cylinder made of an inner core filled with the intracellular space and an outer layer filled with the cell's membrane (Figure 1). Both the intracellular space and the membrane are assumed to be homogeneous such that the study of neuronal electromechanics can be reduced to the study of a simple electromechanical element that we introduce here. Our novel low-dimensional electromechanical model of a neuron couples a spring-dashpot-mass mechanical model of the neuron and an electric circuit model of cell's membrane (Figure 1). Inspired by recent experimental findings by Lu et al. (2006) and Grevesse et al. (2015) we model the macroscopic material neuron as a linear visco-elastic Kelvin–Voigt solid. We use the classic Hodgkin–Huxley equations (Hodgkin and Huxley, 1952) to model the macroscopic electric dynamics of neuron's membrane. The linkage between the Kelvin–Voigt and Hodgkin–Huxley models is achieved by using a Lagrangian formulation and Hamilton's principle as follows.

FIGURE 1

Figure 1. Schematic of the proposed model: the neuron is an axi-symmetric homogeneous circular cylinder whose inner core is the intracellular space (light blue), and the outer layer is the membrane (purple). Due to the symmetry (dashpot line) and material homogeneity, it is enough to study half of the neuron whose properties are encapsulated into a spring-dashpot-mass mechanical system with the spring and dashpot connected in parallel (Kelvin–Voigt model), and the membrane is represented as an electric circuit governed by the classic Hodgkin–Huxley equations.

We start by introducing a Lagrangian of the form:

$\begin{array}{cc}\begin{array}{l}{ℒ}=\frac{1}{2}M{\dot{u}}^2+\frac{1}{2}\tilde{a}{\dot{m}}^2+\frac{1}{2}\tilde{b}{\dot{n}}^2+\frac{1}{2}\tilde{c}{\dot{h}}^2-\frac{1}{2C\left(u\right)}{q}_C{}^2\hfill \\ -\frac{1}{2}k\left(m,n,h\right)u^2,\hfill \end{array}& \left(1\right)\end{array}$

where M is half of the constant mass of the neuron of constant cross-sectional area A, u(t) is the macroscopic (cell level) displacement that depends on time t, k(m,n,h) is the macroscopic spring constant (can be interpreted as a “rescaled” stiffness, as we will show in the results section), C(u) is the macroscopic capacitance of membrane's lipid bilayer modeled as a capacitor of electric charge q_C, and $\tilde{a},\tilde{b},$ and $\tilde{c}$ are positive constants with physical units of Joules. Lastly, m(t),n(t), and h(t) are time-dependent variables between 0 and 1 representing the activations of the Na+ and K+ channels and, respectively, the inactivation of Na+ channel. For simplicity, we denote by $\dot{f}=\frac{df}{dt}$ the first order time derivative of a generic function f(t). The first term of the Lagrangian ${ℒ}$ in Formula (1) is the macro-kinetic mechanical energy, while the second, third and fourth terms are micro-kinetic electric energies. The fifth term in Formula (1) represents a macro-potential electric energy and the last term of ${ℒ}$ is a macro-potential mechanical energy.

Following the variational formulations for electric circuits (Ober-Blöbaum et al., 2013) and for neurons (Chuankui, 2012) we introduce q_Na,q_K, and q_l, the electric charges of Na⁺, K⁺ and leakage channels, respectively. The law of charge conservation provides the following constraint:

$\begin{array}{cc}{q}_C+q_{Na}+q_K+q_l=0.& \left(2\right)\end{array}$

We take q_Na, q_K, q_l, m, n, h, and u as generalized coordinates. By replacing q_C from Formula (2) into the Lagrangian expression (1) we can calculate the variation of the Lagrangian ${ℒ}$ as follows:

$\begin{array}{cc}\begin{array}{l}\delta {ℒ}=li{m}_{\epsilon \to 0}{ℒ}(q_{Na}+\epsilon \delta q_{Na},q_K+\epsilon \delta q_K,q_l+\epsilon \delta q_l,m+\epsilon \delta m,\hfill \\ n+\epsilon \delta n,h+\epsilon \delta h,u+\epsilon \delta u)\hfill \\ =M\dot{u}\delta \dot{u}+\tilde{a}\dot{m}\delta \dot{m}+\tilde{b}\dot{n}\delta \dot{n}+\tilde{c}\dot{h}\delta \dot{h}+\frac{1}{C}{q}_C(\delta q_{Na}\hfill \\ +\delta q_K+\delta q_l)+\frac{1}{2C^2}\frac{dC}{du}{q}_C{}^2\delta u-ku\delta u\hfill \\ -\frac{1}{2}\left(\frac{\partial k}{\partial m}{u}^2\delta m+\frac{\partial k}{\partial n}{u}^2\delta n+\frac{\partial k}{\partial h}{u}^2\delta h\right),\hfill \end{array}& \left(3\right)\end{array}$

where δq_Na, δq_K, δq_l, δm, δn, δh, and δu are variations of the generalized coordinates.

We further define the virtual work done by non-conservative forces as (Ober-Blöbaum et al., 2013; Cusumano et al., 2015):

$\begin{array}{cc}\begin{array}{l}\delta W=-(R_{Na}{\dot{q}}_{Na}\delta q_{Na}+R_K{\dot{q}}_K\delta q_K+R_l{\dot{q}}_l\delta q_l\hfill \\ +\eta \dot{u}\delta u)+(-E_{Na}\delta q_{Na}-E_K\delta q_K-E_l\delta q_l+F_m\delta m\hfill \\ +F_n\delta n+F_h\delta h+f\delta u)\hfill \end{array}& \left(4\right)\end{array}$

In Formula (4) the terms inside the first set of parentheses represent dissipative forces due to the resistors of resistances R_Na,R_K,R_l in the Hodgkin–Huxley model, and due to the linear dashpot in the Kelvin–Voigt model whose damping coefficient is η (can be interpreted as a “rescaled” dynamic viscosity, as we will show in the results section). The second set of parentheses in Formula (4) contains the following generalized forces: the reverse potentials E_Na, E_K, E_l of the Hodgkin–Huxley model, an externally applied mechanical force f, and forces F_m, F_n, F_h which are work conjugates of the gating variables m,n, and respectively h. The choice of signs in Formula (4) guarantees that the virtual work δW is thermodynamically consistent.

We employ now the non-conservative form of Hamilton's principle:

$\begin{array}{cc}{\displaystyle {\int }_{t_1}^{t_2}\left(\delta {ℒ}+\delta W\right)}dt=0,& \left(5\right)\end{array}$

where the variations δq_Na, δq_K, δq_l, δm, δn, δh, and δu are independent and vanish at the arbitrary times t₁, t₂. By replacing Formulas (3) and (4) into the Hamilton's principle, Equation (5), using integration by parts, the independence of the variations δq_Na, δq_K, δq_l, δm, δn, δh, δu and the fact that these variations are zero at t₁,t₂, we obtain the following Euler-Lagrange differential equations:

$\begin{array}{cc}M\ddot{u}+\eta \dot{u}+ku-\frac{1}{2}\frac{dC}{du}{V}^2=f& \left(6\right)\end{array}$

$\begin{array}{cc}{R}_{Na}{\dot{q}}_{Na}=V-E_{Na}& \left(7\right)\end{array}$

$\begin{array}{cc}{R}_K{\dot{q}}_K=V-E_K& \left(8\right)\end{array}$

$\begin{array}{cc}{R}_l{\dot{q}}_l=V-E_l& \left(9\right)\end{array}$

$\begin{array}{cc}\tilde{a}\ddot{m}+\frac{1}{2}\frac{\partial k}{\partial m}{u}^2=F_m& \left(10\right)\end{array}$

$\begin{array}{cc}\tilde{b}\ddot{n}+\frac{1}{2}\frac{\partial k}{\partial n}{u}^2=F_n& \left(11\right)\end{array}$

$\begin{array}{cc}\tilde{c}\ddot{h}+\frac{1}{2}\frac{\partial k}{\partial h}{u}^2=F_h& \left(12\right)\end{array}$

where V = q_C/C is the potential of the capacitor.

Lastly, Kirchhoff's current law needs to be added to the system of Equations (6–12) (Ober-Blöbaum et al., 2013). Replacing Equations (7–9) into Kirchhoff current law yields the well-known Hodgkin–Huxley equation for the membrane potential:

$\begin{array}{cc}C\dot{V}=I-\text{}\text{}\frac{1}{R_{Na}}\left(V-\text{}{E}_{Na}\right)-\text{}\frac{1}{R_K}\left(V-\text{}{E}_K\right)-\frac{1}{R_l}\left(V-\text{}{E}_l\right),& \left(13\right)\end{array}$

where I is a known external current applied on the membrane.

The unknown functions u, V, m,n, and h can be found by solving the coupled Equations (6, 10–13) with appropriate initial conditions. However, in order to solve these equations we need to provide expressions for $F_m,F_n,F_h,\tilde{a},\tilde{b},\tilde{c},C(u),k(m,n,h)$. These expressions are very difficult to prescribe due to insufficient knowledge of neuronal mechanotransduction processes. Thus, for the sake of simplicity, we take f = 0 in Equation (6), and replace Equations (10–12) by the classic Hodgkin–Huxley equations for m, n, h (Dayan and Abbott, 2001):

$\begin{array}{l}\dot{m}={\alpha }_m\left(1-m\right)-{\beta }_{m}m\hfill \\ \dot{n}={\alpha }_n\left(1-n\right)-{\beta }_{n}n\hfill \\ \dot{h}={\alpha }_h\left(1-h\right)-{\beta }_{h}h\hfill \end{array}$

where

$\begin{array}{l}{\alpha }_m=\text{}\frac{0.1(V+40)}{1-exp\left(-0.1(V+40)\right)},{\beta }_m=4exp\left(-0.0556(V+65)\right),\hfill \\ {\alpha }_n=\text{}\frac{0.01(V+55)}{1-exp\left(-0.1(V+55)\right)},{\beta }_n=0.125exp\left(-0.0125(V\text{}+\text{}65)\right),\hfill \\ {\alpha }_h=\text{}0.07exp\left(-0.05(V+65)\right),{\beta }_h=\frac{1}{1+exp\left(-0.1(V+35)\right)}.\hfill \end{array}$

In addition, we take (Dayan and Abbott, 2001):

$\begin{array}{cc}\frac{1}{R_{Na}}=g_{Na}{m}^{3}h\tilde{A},\frac{1}{R_K}=g_K{n}^4\tilde{A},\frac{1}{R_l}=g_l\tilde{A},& \left(14\right)\end{array}$

with à the surface area of the neuron, and g_Na,g_K,g_l the maximal conductances of the Na⁺, K⁺ and respectively leakage currents.

We propose further expressions for C(u) and k(m,n,h). According to Dayan and Abbott (2001) the capacitance is proportional to the surface area of the membrane and since our model is one-dimensional we could for instance assume that the membrane acts like a parallel-plate capacitor. Thus we have:

$\begin{array}{cc}C=c_m\tilde{A}=\frac{ϵ\tilde{A}}{r+u}=\frac{ϵ\tilde{A}}{r\left(1+u/r\right)}\approx \frac{ϵ\tilde{A}}{r}\left(1-\frac{u}{r}\right),& \left(15\right)\end{array}$

where c_m is the specific membrane capacitance, ϵ is membrane's permittivity, and r is the thickness of the membrane. Regarding the expression for the dynamic spring constant k(m,n,h), we hypothesize that the cell stiffens during an action potential. Such an assumption appears to be supported by the observations made by Hille (2001) and Zou et al. (2013). During activation, the neuron stiffens due to the pulling on cytoskeletal elements by the swelling of the cell (Zou et al., 2013), and by the gates in ion channels that act as protein motors (Hille, 2001). According to Formula (14), Na⁺ conductance uses three gates of type m and one gate of type h, while K⁺ conductance uses four gates of type n and we could for instance assume that:

$\begin{array}{cc}k\left(m,n,h\right)=k_0\left(1+m^3(1-h)n^4\right),& \left(16\right)\end{array}$

where k₀ is the spring constant in the inactive state of the neuron.

We observe that in the proposed model the electromechanical coupling is achieved through Equations (6, 10–12), and through Expressions (15–16).

Results

In our simulations we used the following parameters taken from Dayan and Abbott (2001):

$\begin{array}{l}{E}_{Na}=50mV,E_K=-77mV,E_l=-54.387mV,\\ g_{Na}=1.2\frac{mS}{m{m}^2},g_K=0.36\frac{mS}{m{m}^2},g_l=0.003\frac{mS}{m{m}^2}.\end{array}$

The thickness of the membrane is r = 4 nm, the radius of the neuron is r₀ = 2μm (Dayan and Abbott, 2001), an average Young's modulus (stiffness) of the neuron is E₀ = 200Pa (Lu et al., 2006; Zou et al., 2013), half of the neuronal mass is M = 0.1ng (Corbin et al., 2014). The specific membrane capacitance for a neuron in mechanical equilibrium (u = 0) is 0.01 $\frac{\mu F}{m{m}^2}$ and thus from Formula (15) we have:

$c_m=0.01\left(1-\frac{u}{r}\right)\frac{\mu F}{m{m}^2}.$

We also used a value of μ = 4mPa·s for the dynamic viscosity of the neuron (this value was found by Park et al. (2010) for non-neuronal cells).

Under the assumption that the neuron has a circular cylindrical shape, the cross-sectional area is $A=\pi r_0{}^2$. Then the spring constant of an inactive neuron is calculated from equating two different representations of the restoring force in a linear elastic spring: $k_{0}u=E_0\frac{u}{r_0}A$. Thus $k_0=E_0\frac{A}{r_0}$. Similarly, the damping constant is calculated from the shear force to be: $\eta =\mu \frac{A}{r_0}$. Lastly, in all numerical simulations we applied a constant external current per unit (surface) area of 0.1 $\frac{\mu A}{m{m}^2}$.

Because of the numerical stiffness of Equation (6) we solved instead $M\ddot{u}+\eta \dot{u}+ku=0$, with k given by Formula (16), and we solved the classic Hodgkin–Huxley equations with C given by Formula (15). We notice that this simplification preserves a weaker coupling between the mechanical and electrical behaviors of the neuron expressed by Formulas (15) and (16).

We re-wrote $M\ddot{u}+\eta \dot{u}+ku=0$ as a system of first order differential equations:

$\begin{array}{cc}\dot{u}=d,\dot{d}=-\frac{\eta }{M}d-\frac{k}{M}u& \left(17\right)\end{array}$

and used Matlab built-in function ode15s that solves stiff ordinary differential equations. The function ode15s uses a modified linear multistep backward difference formula of order up to five known to have good stability and changes the stepsize of the discretization adaptively according to a numerical scheme that calculates relative and absolute error tolerances (Shampine and Reichelt, 1997).

The Hodgkin–Huxley equations were solved with the following initial conditions:

$\begin{array}{l}V\left(0\right)=-65mV,m\left(0\right)=\frac{{\alpha }_m\left(V(0)\right)}{{\alpha }_m\left(V(0)\right)+{\beta }_m\left(V(0)\right)},\\ n\left(0\right)=\frac{{\alpha }_n(V\left(0\right))}{{\alpha }_n\left(V(0)\right)+{\beta }_n(V(0))},h\left(0\right)=\frac{{\alpha }_h(V(0))}{{\alpha }_h\left(V(0)\right)+{\beta }_h(V(0))}.\hfill \end{array}$

We solved System (17) using two sets of initial conditions:

$\begin{array}{l}\text{Set}1:u\left(0\right)=1nm,d\left(0\right)=\dot{u}(0)=10nm/ms\\ \text{Set}2:u\left(0\right)=0,d\left(0\right)=\dot{u}(0)=1nm/\mu s\end{array}$

Working with Matlab's default values for the relative error tolerance (10⁻³) and the absolute error tolerance (10⁻⁶), the function ode15s solved System (17) and the classic Hodgkin–Huxley equations coupled by Formulas (15) and (16) using a minimum (maximum) stepsize of 0.0051 ms (1.0076ms) for the initial conditions in Set 1, and a minimum (maximum) stepsize of 0.00075 ms(0.0067ms) for the initial conditions in Set 2. In Figures 2, 3 we show the evolutions of the voltage, gaiting variables, displacement and Young's modulus for Set 1 respectively, Set 2. For initial conditions in Set 1, we observe that the action potentials occur and the Young's modulus and the displacement variations appear to be physically admissible and possibly within a healthy range. The dynamics of the stiffness of a neuron is in agreement with the experimental observations in the normally functioning regime reported by Zou et al. (2013). However, for the initial conditions in Set 2 which mimic a more serious traumatic event, not only that the action potentials do not happen anymore (Figure 3A) but also sustained big oscillations of the displacement field are noticed (Figure 3C). In this case the gaiting variables (Figure 3B) as well as the Young's modulus (Figure 3D) remain almost constant. The solutions obtained using the initial conditions in Set 2 might show damaging effects of a very fast initial speed (jabbing) on the material structure and electro-chemical activity of a neuron. As it is apparent from Figures 2, 3, the proposed model is able to capture neuronal mechanotransduction. In Figure 4 we show the evolutions of the displacements obtained using the two sets of initial conditions during 1 ms. While the oscillations are quickly attenuated for Set 1 of initial conditions (Figure 4A) which allows the action potential to develop soon afterwards, for the initial conditions of Set 2 the amplitudes of the oscillations of the displacement field are much higher than in the previous case and do not appear to diminish in time. Also, the membrane's depolarization does not occur in this case (Figure 4B).

FIGURE 2

Figure 2. Results for initial displacement of 1 nm and initial speed of 10 nm/ms: (A) voltage, (B) functions n (blue line), m (red line), and h (green line), (C) displacement, and (D) Young's modulus.

FIGURE 3

Figure 3. Results for zero initial displacement and initial speed of 1 nm/μs: (A) voltage, (B) functions n (blue line), m (red line), and h (green line), (C) displacement, and (D) Young's modulus.

FIGURE 4

Figure 4. A closer look at the displacement evolution during 1 ms in the following cases: (A) initial displacement of 1 nm and initial speed of 10 nm/ms; (B) zero initial displacement and initial speed of 1 nm/μs.

Conclusions and Future Directions

In this paper we proposed a new electromechanical model that couples the mechanical behavior and electro-chemical activity of a neuron and investigated neuronal mechanotranduction through numerical simulations. The neuron was modeled as a liner-viscoelastic Kelvin–Voigt solid whose electro-chemical activity was described by the classic Hodgkin–Huxley equations. We used a Lagrangian formulation and Hamilton's principle to obtain the coupled equations of motion. This approach has the advantage that it can link macroscopic (cell level) as well as microscopic (ionic level) mechanical and electrical information and thus it can describe neuronal mechanotransduction. In addition we assumed that the membrane's capacitance depends on the mechanical displacement of the neuron, while the Young's modulus of the neuron depends on the gating variables in the Hodgkin–Huxley model. Our numerical simulations were done in Matlab using the built-in function ode15s to solve a simplified version of our differential equations. When a constant external electric current was applied and the initial displacement and speed were of orders of magnitude comparable to the size of the neuron, the action potentials occurred and looked similar to the ones observed in healthy neurons. In this case the dynamics of the neuron's stiffness appeared to be in agreement with experimental measurements done on healthy neurons (Zou et al., 2013). However, for very fast initial speeds which could model a serious traumatic event and in the presence of a constant applied external current, high persisting oscillations in the mechanical displacement of the neuron were observed and the action potentials did not happen, suggesting possible structural and functional damage of the neuron.

One of the limitations of the proposed model is coming up with physically relevant expressions that couple capacitance and displacement and respectively stiffness and gating variables, because there are no experimental observations that could guide us. However, given the simplicity of the proposed model, we hope that our approach will inspire future experimental work that will provide empirical relationships among the model's mechanical and electrical parameters. Another limitation of our approach is the use of the Matlab built-in function ode15s to solve the proposed system of stiff differential equations. Shampine and Bogacki (1989) advised caution in drastically reducing the stepsize in the discretization implemented in ode15s since this action may in fact increase numerical error and cause instabilities in the solutions.

In our future work we plan to develop a better numerical solver that will allow us to solve the fully coupled differential equations which are numerically stiff. In addition, we will perform the bifurcation analysis of the model. Lastly, we intend to incorporate in our model ion transport through cell's membrane using the variational formulation from Eisenberg et al. (2010).

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Chuankui, Y. (2012). “A neuron model based on Hamilton principle and energy coding,” in Proceedings of the 2011 2nd International CACS, AISC, Vol. 145 (Berlin; Heidelberg: Springer), 395–401.

Corbin, E. A., Millet, L. J., Keller, K. R., King, W. P., and Bashir, R. (2014). Measuring physical properties of neuronal and glial cells with resonant microsensors. Anal. Chem. 86, 4864–4872. doi: 10.1021/ac5000625

PubMed Abstract | CrossRef Full Text | Google Scholar

Cusumano, J. P., Roy, A., and Li, Q. (2015). Damage dynamics, rate laws, and failure statistics viat Hamilton's principle. Meccanica 50, 77–98. doi: 10.1007/s11012-014-0055-2

CrossRef Full Text | Google Scholar

Dayan, P., and Abbott, L. F. (2001). Theoretical Neuroscience Computational and Mathematical Modeling of Neural Systems. Cambridge, MA: The MIT Press.

Google Scholar

Dickel, G. (1989). Hamilton's principle of least action in nervous excitation. J. Chem. Soc. Faraday Trans. 85, 1463–1468. doi: 10.1039/f19898501463

CrossRef Full Text | Google Scholar

Drapaca, C. S., and Fritz, J. S. (2012). A mechano-electro chemical model of brain neuro-mechanics: application to normal pressure hydrocephalus. Int. J. Numer. Anal. Model. B 3, 82–93.

Google Scholar

Eisenberg, B., Hyon, Y.-K., and Liu, C. (2010). Energy variational analysis of ions in water and channels: field theory for primitive models of complex ionic fluids. J. Chem. Phys. 133, 104104. doi: 10.1063/1.3476262

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldsmith, W. (2001). The state of head injury biomechanics: past, present, and future: part 1. Crit. Rev. Biomed. Eng. 29, 441–600. doi: 10.1615/CritRevBiomedEng.v29.i56.10

PubMed Abstract | CrossRef Full Text | Google Scholar

Goriely, A., Geers, M. G. D., Holzapfel, G. A., Jayamohan, J., Jerusalem, A., Sivaloganathan, S., et al. (2015). Mechanics of the brain: perspectives, challenges, and opportunities. Biomech. Model. Mechanobiol. doi: 10.1007/s10237-015-0662-4. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Grevesse, T., Dabiri, B. E., Parker, K. K., and Gabriele, S. (2015). Opposite rheological properties of neuronal microcompartments predict axonal vulnerability in brain injury. Sci. Rep. 5:9475. doi: 10.1038/srep09475

PubMed Abstract | CrossRef Full Text | Google Scholar

Hille, B. (2001). Ion Channels of Excitable Membranes, 3rd Edn. Sunderland, MA: Sinauer Associates Inc.

Hodgkin, A. L., and Huxley, A. F. (1952). A qualitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764

PubMed Abstract | CrossRef Full Text | Google Scholar

Lanczos, C. (1986). The Variational Principles of Mechanics, 4th Edn. Mineola, NY: Dover Publications Inc.

Google Scholar

Lang, G. E., Vella, D., Waters, S. L., and Goriely, A. (2015). Propagation of damage in brain tissue: coupling the mechanics of oedema and oxygen delivery. Biomech. Model. Mechanobiol. doi: 10.1007/s10237-015-0665-1. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y. B., Franze, K., Seifert, G., Steinhauser, C., Kirchhoff, F., Wolburg, H., et al. (2006). Viscoelastic properties of individual glial cells and neurons in the CNS. Proc. Natl. Acad. Sci. U.S.A. 103, 17759–17764. doi: 10.1073/pnas.0606150103

PubMed Abstract | CrossRef Full Text | Google Scholar

Ober-Blöbaum, S., Tao, M., Cheng, M., Owhadi, H., and Marsden, J. E. (2013). Variational integrators for electric circuits. J. Comput. Phys. 242, 498–530. doi: 10.1016/j.jcp.2013.02.00

CrossRef Full Text | Google Scholar

Paninski, L. (2006). The most likely voltage path and large deviations approximations for integrate-and-fire neurons. J. Comp. Neurosci. 21, 71–87. doi: 10.1007/s10827-006-7200-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, K., Larry, J. M., Kim, N., Li, H., Jin, X., Popescu, G., et al. (2010). Measurement of adherent cell mass and growth. Proc. Natl. Acad. Sci. U.S.A. 107, 20691–20696. doi: 10.1073/pnas.1011365107

PubMed Abstract | CrossRef Full Text | Google Scholar

Shampine, L. F., and Bogacki, P. (1989). The effect of changing the stepsize in linear multistep codes. SIAM J. Sci. Stat. Comput. 10, 1010–1023. doi: 10.1137/0910060

CrossRef Full Text | Google Scholar

Shampine, L. F., and Reichelt, M. W. (1997). The MATLAB ODE suite. SIAM J. Sci. Comput. 18, 1–22. doi: 10.1137/S1064827594276424

CrossRef Full Text | Google Scholar

Wilson, M. T., and Steyn-Ross, D. A. (2008). Subthreshold dynamics of a single neuron from a Hamiltonian perspective. Phys. Rev. E 78:061908. doi: 10.1103/PhysRevE.78.061908

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, S., Chisholm, R., Tauskela, J. S., Mealing, G. A., Johnston, L. J., and Morris, C. E. (2013). Force spectroscopy measurements show that cortical neurons exposed to excitotoxic agonists stiffen before showing evidence of bleb damage. PLoS ONE 8:e73499. doi: 10.1371/journal.pone.0073499

PubMed Abstract | CrossRef Full Text | Google Scholar