Electro-mechanical response of a 3D nerve bundle model to mechanical loads leading to axonal injury
11 Abstract β Objective:
Traumatic brain injuries and damage are major causes of death and disability. We propose a 3D fully coupled electro-mechanical model of a nerve bundle to investigate the electrophysiological impairments due to trauma at the cellular level.
Methods:
The coupling is based on a thermal analogy of the neural electrical activity by using the finite element software Abaqus CAE 6.13-3. The model includes a real-time coupling, modulated threshold for spiking activation and independent alteration of the electrical properties for each 3-layer fibre within a nerve bundle as a function of strain.
Results:
Results of the coupled electro-mechanical model are validated with previously published experimental results of damaged axons. The cases of compression and tension are simulated here, to induce (mild, moderate and severe) damage at the nerve membrane of a nerve bundle, made of four fibres. Changes in strain, stress distribution, and neural activity are investigated for myelinated and unmyelinated nerve fibres, by considering the cases of an intact and of a traumatized nerve membrane.
Conclusion:
A fully coupled electro-mechanical modelling approach is established to provide insights into crucial aspects of neural activity at the cellular level due to traumatic brain injury.
Significance:
One of the key findings is the 3D distribution of residual stresses and strains at the membrane of each fibre due to mechanically-induced electrophysiological impairments, and its impact on signal transmission.
Index Terms β Coupled electro-mechanical modelling, finite element modelling, equivalences, diffuse axonal injury, trauma. I. I NTRODUCTION
Traumatic brain injury (TBI) is caused by mechanical loading to the head due to a sudden acceleration or a blast wave, for example, causing pathologies which range from focal damage of brain tissue to widespread axonal injury [1], [2]. TBI in humans may result from falls, vehicle accidents, sport injuries, military incidents, etc. [1], [3], [4]. The diffuse form of TBI is called Diffuse Axonal Injury (DAI), i.e. a mechanical pathogenesis of an axonal injury, initiated in the deep white matter regions of the brain at the moment of injury [1], [3]. DAI is the most common pathological feature of the mild and severe cases of TBI at cellular and subcellular levels [1], [5], with progressive course [3], responsible for long-lasting neurological impairments following TBI, and high rates of mortality [5]β[7]. Pathological features of DAI include a wide-ranging variety of tissue lesions of the white matter (such as swelling of axonal mitochondria, thinning of the intermodal axoplasm, and development of myelin intrusions [8]), focal haemorrhages, contusions and other brain injuries [1], [3], [7]. Although DAI is classified as an independent category of disease from TBI, the pathological mechanism of DAI is very complex and, currently, there are no diagnostic criteria [3]. A better understanding of
All the authors are with National University of Ireland, Galway (Rep. of Ireland), (email : [email protected]; [email protected]; [email protected] ; [email protected]; )
Electro-Mechanical Response of a 3D Nerve Bundle Model to Mechanical Loads Leading to Axonal Injury
I. Cinelli, M. Destrade, M. Duffy , and P. McHugh mechanically-induced electrophysiological impairments and damage associated with morphological changes of neural cells is urgently needed to improve diagnosis, clinical treatments and prognosis [3], [7]. A number of experimental models of axonal injury reproduce DAI by experimental traumatic insult, where a Traumatic Axonal Injury (TAI) induces DAI to investigate the relationship between mechanical forces and structural and functional responses of axons in experiments [7]. Experimental studies conducted on single axons [9] and nerve fibres [10] aim to induce TAI by applying pressure [11], [12], displacement [1], [10], strain [9], [13], shear strain [14] and electroporation [15]. Although different types of loads seem to initiate TAI, recent studies have shown that the degree of neuronal impairment is directly related to the magnitude and rate of axonal stretch [1], [9], [10], [16]. Beyond a critical threshold [10], axons appear to be vulnerable to stretch-induced changes that induce morphological changes at the microscale level, increasing axolemma permeability [1], [3], [5]. The alteration of axolemma permeability is evidence of mitochondrial damage and neurofilament compaction [7], [8]. Axonal damage is a common manifestation of DAI. Injury-induced axonal damage involves damage of the axonal cytoskeleton, resulting in a loss of membrane integrity and impairment of axoplasmic transport, leading to changes in electrical signal propagation [4], [9], [17]. The alteration of neural activity in a mechanically-injured nerve is called neurotrauma [9], [17], [18]. Although recent experimental studies highlight complex electro-mechanical phenomena at the nerve membrane layer [19], [20], the injury-induced electrophysiological changes of the electro-mechanical activity in neural cells are poorly understood [7], [9]. Quantifying the induced electro-mechanical changes can help to differentiate severity of injury and to understand the alteration in signal propagation at cellular and nerve bundle levels. Computational electro-mechanical models, coupling mechanics and electrophysiology, are powerful tools to investigate and evaluate neurophysiological, neuropathological processes and neurocognitive damage associated with DAI due to injury at the macroscale [1], [2]. This work presents a novel approach for evaluating and quantifying the changes in neural activity due to TAI by using finite element (FE) models. Previous modeling efforts have simulated one-dimensional damage of a nerve fibre [2] and two-dimensional axonal injury of brain tissue [1]. In this work, using the FE software Abaqus CAE 6.13-3, advanced three-dimensional (3D) models explain the link between TBI and DAI at the microscale level [1], [2], by considering the case of TAI at the axonal and bundle levels. Our 3D FE model of a nerve bundle includes a representation of nervous cells made of extracellular media (ECM), membrane, and intracellular media (ICM). As in [21], these components have finite thicknesses, so that the resulting model is a 3-layer nerve bundle. The bundle model here is a section of an idealized geometry of a nerve bundle, which consists of four parallel cylindrical unmyelinated or myelinated fibres, see Fig. 1. In this paper, a series of mechanical loads (such as pressure [11] and displacement [10]) are applied to the bundle to induce a certain level of damage at the nerve membrane of a fibre, altering the fibre activation dynamics and transmission [16], [17]. This model presents a unique framework for investigating how changes in strain and stress distributions alter the function of 3D myelinated and unmyelinated nerve fibres and bundles. In contrast to previous studies [1], [2], this model includes full coupling between mechanical and electrical domains [22], [23] where the neural electrical activity, as described in [24], is coupled to the mechanical domain through piezoelectricity [20] and electrostriction [19], reproducing biophysical phenomena accompanying the action potentials as observed experimentally [20], [25]. Here, the electrical activity is directly coupled to mechanical deformation by using electro-thermal equivalences in a coupled thermo-mechanical FE analysis [26]. Electro-thermal equivalences and equivalent material properties have been shown to provide an efficient approach to resolve 3D electrical problems in a coupled electro-mechanical analysis in Abaqus CAE 6.13-3 [26]. This paper builds on the work reported in [22], [23]. Here, we report details of how electro-mechanical coupling is represented through thermo-mechanical equivalence; the nerve bundle model and its geometry; and the damage criterion and its implementation. The electro-mechanical coupling approach is validated based on the strain-based damage criterion [2], [16], [22]. This criterion refers to traumatic mechanically-induced damage on a nerve fibre (refer to [2], [16]) in which the resting ionic potentials of the Hodgkin and Huxley (HH) model are shifted by , simulating experimental evidence of damaged nerve fibres [16], [17], as in [2]. The injury threshold takes into account axonal strain along the nerve fibre length only [1], [2]. The strain along the fibre length has been shown to be a physiologically relevant injury criterion for multiscale TBI models [1], [2]. Two cases, (i) dynamic pressure and (ii) displacement loads at the bundle level, in which only one fibre is activated, are considered. This approach has the potential to generate useful insights in studying the mechanics behind neuro-physiology, as observed experimentally in damaged nerve membranes of clinical cases (such as multiple sclerosis) [9], [18], [27].
Fig. 1. (a): frontal view and (b): isometric view of the 3-layer nerve bundle, made of four fibres. Fibre and the Ranvierβs node length is , while the radial thickness of the layer is equal to [26], [28].
II. M ETHOD A. Electro-mechanical Coupling
The electro-mechanical coupling of the action potential [24] is achieved through modelling of the nerve membrane as a piezoelectric material [19], building on the equivalent electro-thermal modelling approach for predicting the electrical behaviour of nervous cells described in [21], [22], [26]. This can be described in terms of the piezo-elastic relation given in equation (1) where πΊ is the total strain vector, π· is the compliance matrix, π is the mechanical stress vector, πΉ is the piezoelectric strain coefficient vector, β is the thickness of the piezoelectric layer, and βπ is the voltage difference across it [29], [30]. πΊ = π·π + πΉ(βπ/β) (1) Similarly, the thermo-elastic strain-stress relation is described in (2), where πΆ is the thermal expansion coefficients vector and π₯π is the temperature difference [29], [30]. πΊ = π·π + πΆβπ (2) By the electro-thermal analogy [22], [26], [29], [30], the electric field, approximated here as the voltage across the membrane divided by its thickness, is equivalent to a thermal load, while the piezoelectric constants are equivalent to the thermal expansion coefficients, see (1) and (2). Simulating the nerve membrane as the dielectric component of a parallel plate capacitor [24], [29], the displacement follows the gradient of the electric field using the approach presented in [29], [30], see (2). That is, the piezoelectric effect is only relevant in the through-thickness direction, represented here with orthotropic piezoelectric constants, approximately per [20] in the thickness direction and zero in the longitudinal and circumferential directions. Meanwhile, incompressible isotropic mechanical properties [25] are assumed. As the electrical depolarization is initiated in the nervous cell [24], compressive forces on the nerve membrane arise from piezoelectricity [20], see (2), leading to changes in shape (electrostriction [20], [25], [31]) and corresponding changes in capacitance. The electrical capacitance per unit area, πΆ , changes as the square of the voltage [25], [31], see (3) and Fig. 2. Here, π is the fractional increase in capacitance per square volt (equal to !" [31]) and βπ is the surface potential difference (about β70 ππ [31]). πΆ(π) = πΆ(0)[1 + π(π + βπ) " ] (3) Hence, this approach establishes a fully electro-mechanical coupling accounting, at the same time, for the effect of mechanical loads on electrical response and the effect of the electrical activity on the nerve structure. B. Model
The bundle model simulates the exchange of charges in four identical fibres as shown in Fig. 1. Each fibre consists of a cylindrical region of ICM enclosed by a thin membrane and surrounded by a region of ECM [26]. Two bundle models are used in this study: a fully unmyelinated bundle and a fully myelinated bundle. The neurite radii are: π = 0.477 Β΅π , π % = 0.480 Β΅π and π &$% = 0.500 Β΅π [26] for the ICM, membrane and ECM layers, respectively. As a first step, this analysis is focused on the radial distribution of charges rather than on longitudinal variations. The length of the bundle is for a diameter equal to in both cases, see Fig. 1, within the range of the human optic axon [32]. To model a myelinated fibre, we assume an ICM enclosed by a single layer, which is periodically-partitioned along the fibre length, similarly to the histologic section of a myelinated fibre, see Fig. 1 (d.1) and (d.2). The insulation sheath of myelin around the nerve fibre is modelled as an insulating layer, which replaces the membrane layer at regular intervals along the fibre [28], see Fig. 1 (d.1). Different conductivity values are assigned to denote the myelin and membrane sections [28], see Fig. 1 (d.2). The width of the piecewise conductive membrane regions (or Ranvierβs nodes) is and the internode distance is [28], see Fig. 1 (d.2). Incompressible isotropic mechanical properties [25] are assumed in both models, with Youngβs Modulus equal to and a Poissonβs ratio of [25] (close to incompressibility). Then, the electrical model parameters for unmyelinated and myelinated fibres are taken from [33] and [2], respectively. The resistance per unit length of ICM and ECM is , and for the myelin it is . Then, the capacitance per unit area is " β for the ICM, " β for the ECM, and !β ππΉ ππ " β for the myelin layer. The corresponding value of the membrane are the HH model values [24] for the steady state and are dependent on voltage, strain and spatial coordinates during loading conditions, see section D. According to previous published work [34], [35], a neurite refers to a part of the body of a nervous cell (axon) that consists of two structurally distinct regions (axoplasm and membrane). In contrast, a fibre usually refers to an axon (as axoplasm and membrane) with a myelin layer. In this work, neurite refers to myelinated fibre too, because the myelin layer is modelled as periodical insulating layers at the nerve membrane layer, whose radial extension has been neglected, as in [28]. C. Damage Evaluation
Strain-based damage affects sodium and potassium reversal potentials (πΈ () (π * ) and πΈ + (π * )) , see (4), of the HH model, simulating the changes in ionic concentration across the nerve membrane depending on the membrane strain (π * ) [2]. The resting potentials at physiological conditions are πΈ (), and πΈ +, [2], [24] see (4). If a maximum strain, πΜ , is exceeded, then the reversal potentials are zero, otherwise the changes follow the equation (4) if π * < πΜ [2]: SπΈ () (π * ) = πΈ (), (1 β (π * πΜβ ) - )πΈ + (π * ) = πΈ +, (1 β (π * πΜβ ) - ) (4) Here, the strain threshold, πΜ , is set at as an indicator of the onset of functional damage at which there is a probability of having morphological injury during an uniaxial displacement test of a nerve fibre [10]. The parameter πΎ , equal to [2], is an index referring to the sensitivity of the damage to small versus large deformation, see [2]. Additionally, the reversal potential of the leak ions πΈ . ! is not influenced by the strain but varies based on changes in gradient concentrations of potassium and sodium across the membrane [2], see Fig. 2. The changes for the leak ions can be derived from the sum of the total membrane currents at resting conditions [2], [24]. D. Implementation
The implementation of the coupled Hodgkin and Huxley (HH) model is shown in Fig. 2 (on the right), in contrast to the uncoupled HH model (shown on the left). By using the electro-thermal equivalence [22], [26], the implementation of the neural activity, the distribution of voltage, and the generated strain can be seen in 3D by using well established coupled thermo-mechanical software simulation tools. In the coupled model [22], the membrane electrical conductivity changes in response to the membrane voltage produced during the action potential, incorporating the effects of strain [2], [24], according to (4), while the electrical capacitance per unit area changes with the square of the voltage [31], see (3). The HH reversal voltage potentials change due to the strain at the membrane [2], [16], hence the threshold of spike initiation changes as in [24]. In all cases, strain is caused by membrane piezoelectricity as described in (1), and external mechanical loading when applied. The model is implemented as a coupled thermo-mechanical model in Abaqus CAE by using user-defined subroutines (DISP, USDFLD and UMATHT) to assign thermal equivalent electrical properties to the membrane of each fibre, independently, based on the spike initiation [36], strain [2], [16] and voltage [31] generated at the same membrane, see (2) and (3). Electrical model parameters for other regions of the unmyelinated and myelinated fibres are taken from [33] and [2], respectively. E. Boundary Conditions
In all cases, a voltage Gaussian distribution (with zero-mean and 0.4 variance) is the upper-threshold stimulation applied on Fibre or of the transmembrane voltage) assuming an intact or traumatized membrane, respectively, as in [2], [16]. An encastrΓ© boundary condition is enforced at the origin of each model, while no rotation nor shear is allowed. Then, two cases of mechanical loads applied at the bundle have been considered. As a first step to assess this novel coupling, only frequency-independent loading conditions are considered throughout, after the initial steady-state, so the solutions are quasi-static [21]. In the first case of applied mechanical loading, an instantaneous uniform compression is applied to the bundle to simulate injury conditions. Three values of pressure are modelled simulating mild ( less than , moderate (55 β 95 πππ) and severe ( higher than pressure [11]. In the second case, axial strain conditions that reproduce the uniaxial test in [10] are applied. Two values of instantaneous uniform stretch are applied as a displacement boundary condition to simulate and of total axial deformation, π , at which the probability of inducing morphological injury during the elongation test is 5% and 25% respectively [10]. Fig. 2. Flowchart of the code describing the active behaviour of the nerveβs membrane: on the left, the HH dynamics [24] and on the right, the fully coupled HH dynamics. Here, a Gaussian voltage distribution elicits the action potential in a 3D model of a nervous cell. By using electro-thermal equivalences, the HH model is implemented as an equivalent thermal process, in which the membraneβs conductivity changes as in [24] and the capacitance, πΆ ! , changes as in [31]. The HH parameters are changing based on the temperature (here the equivalent of voltage [21]) and strain at the membrane [2]. The strain π generated in the model is a function of temperature, π , and thermal expansion coefficients, see (3) and (4).Voltage, current, strain and stresses distribution are only a few of the 3D results released by Abaqus by equivalence. F. Validation
The traumatic mechanically-induced damage [2], [16] is validated by using a coupled left-shift ( πΏπ ) version of the HH model [16] reproducing experimental evidence [16], [17], as in [2]. The LS denomination refers to the induced shift of the transmembrane potential to positive values, simulating a nerve membrane subjected to mechanical trauma [2]. It is called coupled LS to highlight the close link of mechanical strain on the electrical properties of the nerve membrane. In contrast to [24], the ionic currents of the modified HH model include the fraction of affected ionic channels by the strain ( π΄πΆ ) and coupled left-shift ( πΏπ ) variables, see (5) and (6) [37]. When a trauma occurs at the nerve membrane, the ionic resting potentials undergo a voltage-shift ( πΏπ ), here, equal to , as a reference value for trauma-induced kinetic changes observed in experiments [37]. This allows to simulate uniaxial loading experiments [13], [38] of a nerve bundle, by implementing a pre-imposed shift in the HH dynamics, simulating the impact of applied mechanical strain at the nerve membrane. Then, the integrity of the nerve membrane is modelled as the fraction of the voltage-time dependent activation and inactivation variables (i.e. π, π and β [24]) of the HH model, shifted by πΏπ (here, indicated as π /0 , π /0 and β /0 ) , see (5) and (6) [37]. The model accounts for the case where only a fraction of nodal channels ( π΄πΆ ) undergoes a πΏπ , while the rest of the membraneβs channels, ( , remain intact [16]. Here, only the extreme cases of the entire membrane being traumatized (π΄πΆ = 1) or intact (π΄πΆ = 0) are shown as illustrative examples. In (5), the sodium current, πΌ () , and, in (6), the potassium current, πΌ + , [16]. The membrane potential is π * , and the sodium, πΜ () , and potassium, πΜ + ,conductances are, respectively, equal to " _ and
36 ππ ππ " _ [24]. This model is validated considering two cases of interest and two approaches. The first case (I) is the nerve membrane at physiological conditions, i.e. intact ( πΏπ = 0 ) and non-traumatized (
π΄πΆ = 0 ). In second case (case II), the membrane loses integrity (
π΄πΆ = 1) due to damage (πΏπ = 20 ππ).
Then, in the first approach, the reversal ionic potentials in (5)-(6), πΈ () and πΈ + , are constant values, equal to and β77 ππ respectively [24], as in [16]. Here, the signals are changing because of the changes in conductance produced by changes in the voltage-dependent variables (π /0 , π /0 and β /0 ) . In a second approach, the reversal potentials are changing with πΏπ , see (4), as in [2]. Here, the signals are changing due to the changes in both the conductance and the reversal ionic potentials. Fig. 3 shows the changes in action potential due to the left-shift voltage implemented as in [16] (indicated with *) and as in [2] (indicated as **). Those two conditions are the ones that can be directly compared to [2], [16]. All the equations are collected in the Appendix in greater detail. III. R ESULTS A. Model Analysis
The electrophysiological impairments associated with structural damage of the neuron mechanics affect the electrostriction and the electric field across the membrane, (3), while the piezoelectric properties of the nerve membrane are assumed to be constant, (2), [20], [25]. The results in Fig. 3 confirm the ability of the model to reproduce the membrane voltage and ionic current waveforms of a damage-induced injury. Here, the coupled LS model of HH reproduces the changes in hyperpolarization in a damaged axon, left-shifting the affected ionic channels by , as in slow-severe cases, as in [2], [16]. Fig. 3 refers to Fibre πΌ . ! , potassium, πΌ + , and sodium, πΌ () ) vs. time in Cases I, II (*) and II (**). As shown, the mechanically-induced voltage-shift of the resting potentials leads to a time-delay and a lower amplitude of the voltage signal, due to leaky channels, similar to trends observed in [16] and [2]. In Case II (*), the action potential is shifted by and its peak at has a magnitude of about . In Case II (**), the action potential is shifted by and the maximum magnitude is about β46.5 ππ at . To illustrate the utility of the model in predicting the nerve mechanical response, Figs. 4 (a) and 4 (b) show the corresponding displacement vs. time of an unmyelinated nerve membrane and myelinated nerve membrane, respectively, in a damaged bundle along the radial direction (i.e. the π₯ β ππ₯ππ ). Data are taken at πΌ () = [π β β(1 β π΄πΆ) + π /0β β /0 π΄πΆ](π * β πΈ () )πΜ () (5) πΌ + = [π (1 β π΄πΆ) + π /0β π΄πΆ](π * β πΈ + )πΜ + (6) the position of maximum radial displacement on Fibre πΏπ =20 ππ and
π΄πΆ = 1 . In both cases, data are plotted at the time instant when the mechanical displacement is at its maximum.
Fig. 3. (a) the membraneβs potential in Case I, action potential (πΏπ = 0ππ; π΄πΆ = 0) ; Case II, damaged traumatized membrane (πΏπ = 20ππ; π΄πΆ = 1) . In *, the reversal ionic potentials have constant values as in [16], and, in **, they change according to equation (4), as in [2]. In (b), the Current Density [π΄ ππ " β ] on Fibre For each case in Fig. 3, residual compressive forces are found in unmyelinated and myelinated fibres, see Fig. 4 and 5, due to the biophysical activity of the nerve membrane [20], [31]. In an unmyelinated bundle, the displacement peak is β3.3 ππ for Case I, β2.5 ππ for Case II (*), and β0.6 ππ for Cases II (**). In a myelinated bundle, the peak is β0.81 ππ for Case II (*) and β0.27 ππ for Case II (**). The peaks of the radial displacements occur at the membrane potential peaks for each case. In a myelinated bundle, lower displacements are found at the membrane than in the unmyelinated case, see Fig. 4 (b) and 5 (e)-(h). The myelin layer constrains the deformation of the nerve membrane, reducing the displacement by about an order of magnitude on the same unmyelinated fibre, see Fig. 4 (b) and Fig. 5 (h). Assuming the Ranvier node regions are aligned in a bundle, the total displacement of the bundle is driven by the same deformation as the activated fibre, see Fig. 5 (h). B. Mechanical Loading Cases of Interest
In this section, the fully coupled electro-mechanical model of Fig. 2 is applied to investigate different mechanical loading conditions directly using equations (1-4), instead of as a simulated voltage-shift. Pressure Loads
Fig. 6 (a) and Fig. 6 (b) show the hyperpolarization and current densities of an unmyelinated bundle and a myelinated bundle under mild ( ), moderate ( ) and severe ( ) pressures inducing Traumatic Axonal Injuries (TAI) [11]. The case of extreme pressure ( ) is also considered. In contrast to the reference case of an intact nervous cell (
π = 0 πππ and
π΄πΆ = 0 ), all pressure loads reduce the magnitude of the action potential and shift it over time, correlating with results in Fig. 3. Here, the strain applied at the nerve membrane by compressing the bundle shifts the ionic resting potentials of the fully coupled HH model by a quantity which varies depending on the magnitude of the applied load, see (4). Thus, only the π΄πΆ variable is, here, considered. Compression levels in the range of to have a similar impact on the signal transmission both in terms of reduced magnitude and shift over time; this is found to be due to similar strain values read at the nerve membrane, because of the linear elastic assumption. In these cases, the peak of the action potential is about β17 ππ with of time-delay. Here, only slight differences are found for a traumatized vs. non-traumatized nerve membrane (π΄πΆ = 1 vs. 0) when mild-to-severe pressures are applied, and therefore only results for AC = 0 are included, see Fig. 6. Fig. 4. (a): mechanical displacement of the unmyelinated nerve membrane; (b): displacement of a myelinated nerve membrane of the Fibre π₯ β ππ₯ππ ), π . Fig. 5. Frontal and isometric views of the total displacement in (a) an unmyelinated and (e) a myelinated nerve bundle model. In (b) and (f), the ECM; in (c) and (g), the ICM of the two bundle types. In (d), the isometric view of the nerve membrane and, in (h), the Ranvier nodes and the myelin sheath of the myelinated bundle. (a)-(d) and (e)-(h) are Case II (**) [2] applied to an unmyelinated and myelinated nerve bundle. Data are taken at the peak of the action potential in both cases.
In contrast, application of extreme pressure changes the resting potential of ions reducing the current flow across the nerve membrane, see Fig. 6 (b), and reducing the magnitude of the action potential to zero, see Fig. 6 (a). An extreme pressure leads to both a reduction in magnitude and an increase of the voltage baseline up to β24 ππ for an intact membrane, and up to β7ππ for a traumatized membrane, see Fig. 6 (a) and (b). Here, the homeostatic balance of charges across the nerve membrane vanishes with an extreme pressure applied over a traumatized nerve membrane (π΄πΆ = 1) , because the strain levels at the nerve membrane are close to the threshold value assumed in (4). With mild-to-severe pressure loads from to , the contraction of the nerve membrane due to the electrostriction has a similar trend to the case of πΏπ = 0 ππ and
π΄πΆ = 1, see Fig. 3. Figs. 7 (a) and 7 (b) show the radial mechanical displacement on Fibre for the extreme pressure case (not shown here) and by for the mild-to-severe cases [11], respectively, see Figs. 7 (a) and 7 (b). Higher deformation at the nerve membrane changes the ionic resting voltages leading to higher mechanical displacements, see (4) and Fig. 8. Fig. 6. (a) Membrane Potential [ ππ ] and (b) Current Density [π΄ ππ " β ] on Fibre ), moderate ( ) and severe ( ) pressures [11]. The extreme case of π = 1πΊππ is also considered. AC is the fraction of affected ionic channels by the strain:
π΄πΆ = 0 is for an intact membrane and
π΄πΆ = 1 for a traumatized membrane [16]. Data are the maximum radial displacement of a node on Fibre
Then, for the mild to severe cases of TAI-induced pressure [11], the resting voltage potentials are changed due to the induced deformation in the bundle and the magnitude of the action potential is, hence, reduced [2], see Fig. 7. The peak of the action potential is higher in a compressed unmyelinated bundle than in a compressed myelinated one. Within the range of pressure levels considered, an unmyelinated layer displaces according to the charges exchanged across the nerve membrane, see Fig. 7. On Fibre β1 ππ , β1.10 ππ , and β1.40 ππ , respectively, see Fig. 7 (a). In a myelinated bundle, the charge-induced displacement of a myelinated fibre is much less than in an unmyelinated bundle and therefore its displacement is more in response to the loading condition than to electrostriction, see Fig. 7 (b). Here, on Fibre β0.14 ππ , β0.40 ππ , and β1.20 ππ respectively, see Fig. 7 (b). This model shows greater membrane displacements in an unmyelinated fibre than in a myelinated fibre, see Fig. 7 (a) and 7 (b). The myelin layer constrains the deformation of the Ranvier nodes, which are the only regions throughout the fibre to show voltage-induce membrane displacement [20]. At the nodes, the applied compression acts in opposition to the electrostriction, because of the negative value of the action potential. Fig. 7. Radial displacement [ππ] of (a) an unmyelinated bundle and (b) a myelinated bundle. Uniform applied pressures are classified as mild ( ), moderate ( ) and severe ( ) pressures [11]. Data are the maximum radial displacement of a node on Fibre Displacement Loads
In Figs. 8 (a)-(b) and 8 (c)-(d) a displacement boundary condition applied along the length of the unmyelinated and myelinated bundles, respectively, simulates and of total deformation [7]. In an unmyelinated bundle, see Fig. 8 (a)-(b) and Fig. 9 (a)-(d), the action potential is shifted through time durations by about for π = 5% and π΄πΆ = 0 , and by for π = 5% and
π΄πΆ = 1 . For higher deformation, the action potential is delayed by , showing similar results both with
π΄πΆ = 0 and
π΄πΆ = 1 , representing the loss of ionic gradient across the nerve membrane, see Fig. 8 (b). In contrast, in a myelinated bundle, see Fig. 8 (b)-(d) and 9 (e)-(h), no significant differences have been found between intact and traumatized membranes. Here, for and deformation, the membrane potential is shifted up to β60ππ and β43ππ , respectively, while the peak of the action potential is shifted at and , respectively, see Fig. 8 (c). Lower current density at the Ranvier nodes is mainly due to lower voltage gradient and higher localized strain, see Figs. 8 (d) and 9 (h). IV. D ISCUSSION
In contrast to previous studies [1], [2], this paper shows the advantages of a fully coupled electro-mechanical 3D framework to investigate the details of neural activity, combining real-time fully coupled electro-mechanical phenomena, modulated threshold for spiking activation, and independent alteration of the electrical properties for each fibre in the 3-layer nerve bundle, made of membrane (or Ranvierβs nodes and myelin sheath), ICM and ECM. The electro-mechanical coupling, based on electro-thermal equivalences [22], [23], [26], allows for reliable simulation of changes in electrostriction and neural activity due to mechanical damage, as seen in experiments [10], [11], [16].
Fig. 8. (a)-(b) Membrane Potential [ ππ ] and (c)-(d) Current Density [π΄ ππ " β ] on Fibre and of total deformation π applied [10]. AC is the fraction of affected ionic channels by the strain: π΄πΆ = 0 is for an intact membrane and
π΄πΆ = 1 for a traumatized membrane [16]. Data are taken at the maximum displacement along the bundle middle axis, i.e. π§ β ππ₯ππ , on Fibre Fig. 9. Voltage distribution (ππ11) in an unmyelinated nerve bundle, in (a)-(d), and in a myelinated nerve bundle, in (e)-(h), for elongation. Frontal and isometric view of an unmyelinated and myelinated nerve bundle model in (a) and (b), respectively. In (b) and (f), the ECM, in (c) and (g) the ICM of the two bundle types. In (d), the isometric view of the nerve membrane and, in (h), the Ranvierβs nodes and the myelin sheath of the myelinated bundle. In this study, two cases of interest provide insights into the electrophysiological impairments of axonal injury due to sudden TAI-induced pressures and displacements. Differences in signal transition arise in the bundle for each fibre, depending on the fibre type. In the bundle, Fibre and elongation greater than , see Figs. 7 and 9. Unmyelinated nerve fibres show the greatest changes with mechanical loads due to the higher current density per area exchanged at the membrane. In contrast, lower density of ions per area can flow thought the membrane at the Ranvier nodes, while the myelin layer constrains the electrostriction accompanying the neural activity. Results show that in the myelinated bundle the activation of the fibres is compromised because the fibre is not able to follow the pattern of activation under pressure and displacement loads. These two loading conditions are cases of clinical interest. A constant pressure on an axon is representative of injuries [1], [3] and plaques in demyelinating diseases, acting as conduction block for action potentials of neighbouring axons inhibiting the myelinisation process. Instead, elongation tests might explain the relatively long time period needed for spontaneous membrane repair after fast strain [13]. Finally, as highlighted earlier, the results refer to a cylindrical bundle made of four identical fibres with characteristics within the range of the human optic axon [32]. Although real nerve bundles are made of a higher number of fibres with different calibre [32], the use of a simplified geometry was needed to assess electro-mechanical equivalences in a 3D FE model and to limit the computational cost. On-going work is focused on extending this technique to nerve bundle with different calibre (made of unmyelinated, myelinated and mixed fibres) and multiple activation of fibres. V. C ONCLUSION
We propose a fully coupled electro-mechanical framework for modelling the biophysical phenomena accompanying neural activity. Here, the coupling is based on an electro-thermal analogy by modelling the piezoelectric effect as a thermal expansion phenomenon in Abaqus CAE 6.13-3. This approach allows us to model and generate insights into aspects of neural activity, such as electrostriction and piezoelectricity, and to correlate these with experimental observations. The model, built on previously published work [22], [23], [26], generates a fully 3D simulation of ion channel leaking for nerve fibres under pressure and displacement loads. In conclusion, in this model: β’ Time-shift, signal magnitude and nerve membrane potential baseline are dependent on the total strain, voltage and size of the fibre; β’ Lower strain and lower electrophysiological changes are found in myelinated fibre Vs. unmyelinated fibre; β’ Unmyelinated fibre are needed to share information across the fibre, rather than throughout its length. This model can contribute to the understanding the causes and consequences of TBI and DAI to improve diagnosis, clinical treatments and prognosis by simulating the mechanical changes accompanying the changes in signal transmission in TAI-induced loading conditions. VI. A CKNOWLEDGMENT
The authors gratefully acknowledge funding from the Galway University Foundation, the Biomechanics Research Centre, and the Power Electronics Research Centre, College of Engineering and Informatics, NUI Galway, Ireland. A
PPENDIX
The fully coupled Hodgkin and Huxley model reported in this work is based on the electro-thermal equivalences [39] and on the electro-mechanical coupling, see (A.3) and (A.4). First, (A.1) and (A.2) show the implementation of the 3D Heat Diffusion Equation as Cable Equation, see [39]. The electrical material properties are substituted to the corresponding thermal properties, as discussed in [39]. In the Heat Equation (A.1), π is the temperature [ Β° πΆ or πΎ] , π is the mass density [ππ/π β ] , π is the specific heat capacity [π½/(ππ πΎ)] , π is the thermal conductivity [π/(ππΎ)], and Q is the heat source density [π/π β ]; correspondly, in the Cable Equation (A.2) applied to the nerve, π * is the transmembrane potential [π] , π [1/π] denotes the surface-volume ratio, πΆ * [πΉ/π " ] is the nerve membrane capacitance, π [π/π] is the electrical conductivity of the membrane, and πΌ [π΄/π] is the ionic current . ππ ππ ππ‘β β π» β (ππ»π) + π = 0 (A.1) πΆ * π ππ * ππ‘β β π» β (π π»π * ) + π πΌ = 0 (A.2) Second, (A.3) and (A.4) display the piezo-elastic strain-stress relationship vs. the thermo-elastic strain-stress relationship, respectively (see Method β Section A). In the electro-mechanical coupling, the piezoelectric strains are simulated by using the analogue quantity corresponding to the piezoelectric coefficients, i.e. the expansion coefficients πΆ . πΊ = π·π + πΉ(βπ/β) (A.3) πΊ = π·π + πΆβπ (A.4) Then, simulating quasi-static solutions [39], (A.5) shows the balance of ionic currents modelling the nerve membrane, as in the Hodgkin and Huxley model [24]. πΆ % ππππ‘ = πΜ () π β β(π β πΈ () ) + πΜ + π (π β πΈ + ) + πΜ . (π β πΈ . ! )= πΊ () (π β πΈ () ) + πΊ πΎ (π β πΈ + ) + πΊ . ! (π β πΈ . ! ) = πΌ () + πΌ + + πΌ . ! (A.5) Following trauma, electrical alterations can be due to: (i) the loss of integrity of the membrane [16] where a fraction of affected ionic channels ( π΄πΆ ) is assumed to be traumatized; and (ii) transmembrane voltage shift ( πΏπ ) to positive values [16], representing a different osmotic gradient due to the trauma-induced strains. Sodium, πΌ () , potassium, πΌ + , and leak ions current, πΌ . ! , are shown in (A.6), (A.7) and (A.8), respectively (see Method β Section F). πΌ () = [π β β(1 β π΄πΆ) + π /0β β /0 π΄πΆ](π * β πΈ () )πΜ () (A.6) πΌ + = [π (1 β π΄πΆ) + π /0β π΄πΆ](π * β πΈ + )πΜ + (A.7) πΌ . ! = πΜ . ! (π * β πΈ . ! ) (A.8) Additionally, the reversal potentials of potassium, (A.9), sodium, (A.10), and leak ions, (A.11), as the axonal component of the total strain applied at the bundle is changing (see Method β Section C). The membrane resting potential is π . πΈ () (π * ) = πΈ (), (1 β (π * πΜβ ) - ) (A.9) πΈ + (π * ) = πΈ +, (1 β (π * πΜβ ) - ) (A.10) πΈ . ! = u1 + πΊ () + πΊ + πΊ . ! v π β πΊ () πΈ () (π * ) + πΊ + πΈ + (π * )πΊ . ! (A.11) Finally, the electrical capacitance per unit area is changing with voltage as in (A.12) (see Method β Section A). πΆ(π) = πΆ(0)[1 + π(π + βπ) " ] A.12 R
EFERENCES
1. R. M. Wright and K. T. Ramesh, An axonal strain injury criterion for traumatic brain injury,
Biomech. Model. Mechanobiol. , vol. 11, no. 1β2, pp. 245β260, 2012. 2. A. JΓ©rusalem, J. A. GarcΓa-Grajales, A. MerchΓ‘n-PΓ©rez, and J. M. PeΓ±a, A computational model coupling mechanics and electrophysiology in spinal cord injury,
Biomech. Model. Mechanobiol. , vol. 13, no. 4, pp. 883β896, 2014. 3. J. Ma, K. Zhang, Z. Wang, and G. Chen, Progress of Research on Diffuse Axonal Injury after Traumatic Brain Injury,
Neural Plast. , vol. 2016, p. 7, 2016. 4. Y. P. Zhang, J. Cai, L. B. E. Shields, N. Liu, X. M. Xu, and C. B. Shields, Traumatic Brain Injury Using Mouse Models,
Transl. Stroke Res. , pp. 1β18, 2014. 5. D. H. Smith and D. F. Meaney, Axonal damage in traumatic brain injury,
Porg Clin Neurosci , vol. 6, no. 6, pp. 483β492, 2000.
6. H. C. Wang and Y. Bin Ma, Experimental models of traumatic axonal injury,
J. Clin. Neurosci. , vol. 17, no. 2, pp. 157β162, 2010. 7. L. E. Delisi, Handbook of Neurochemistry and Molecular Neurobiology,
Handb. Neurochem. Mol. Neurobiol. Schizophr. , vol. 27, pp. 107β241, 2009. 8. M. W.L. and G. D.I., Loss of Axonal Microtubules and Neurofilaments after Stretch-Injury to Guinea Pig Optic Nerve Fibers,
J. Neurotrauma , vol. 14, no. 9, pp. 603β615, 1997. 9. J. a Galbraith, L. E. Thibault, and D. R. Matteson, Mechanical and electrical responses of the squid giant axon to simple elongation.,
J. Biomech. Eng. , vol. 115, no. 1, pp. 13β22, 1993. 10. A. C. Bain and D. F. Meaney, Tissue-level thresholds for axonal damage in an experimental model of central nervous system white matter injury.,
J. Biomech. Eng. , vol. 122, no. 6, pp. 615β622, 2000. 11. S. Hosmane et al. , Valve-based microfluidic compression platform: single axon injury and regrowth,
Lab Chip , vol. 11, no. 22, p. 3888, 2011. 12. P. E. Gallant, The direct effects of graded axonal compression on axoplasm and fast axoplasmic transport.,
J. Neuropathol. Exp. Neurol. , vol. 51, no. 2, pp. 220β30, 1992. 13. R. Shi and J. Whitebone, Conduction deficits and membrane disruption of spinal cord axons as a function of magnitude and rate of strain.,
J. Neurophysiol. , vol. 95, no. 6, pp. 3384β90, 2006. 14. M. C. LaPlaca, D. K. Cullen, J. J. McLoughlin, and R. S. Cargill, High rate shear strain of three-dimensional neural cell cultures: A new in vitro traumatic brain injury model,
J. Biomech. , vol. 38, no. 5, pp. 1093β1105, 2005. 15. P. E. Gallant and J. A. Galbraith, Axonal structure and function after axolemmal leakage in the squid giant axon.,
J Neurotrauma , vol. 14, no. 11, pp. 811β822, 1997. 16. P. A. Boucher, B. JoΓ³s, and C. E. Morris, Coupled left-shift of Nav channels: Modeling the Na+-loading and dysfunctional excitability of damaged axons,
J. Comput. Neurosci. , vol. 33, no. 2, pp. 301β319, 2012. 17. N. Yu, C. E. Morris, B. JoΓ³s, and A. Longtin, Spontaneous Excitation Patterns Computed for Axons with Injury-like Impairments of Sodium Channels and Na/K Pumps,
PLoS Comput. Biol. , vol. 8, no. 9, 2012. 18. D. M. Geddes, R. S. Cargill, and M. C. LaPlaca, Mechanical stretch to neurons results in a strain rate and magnitude-dependent increase in plasma membrane permeability.,
J. Neurotrauma , vol. 20, no. 10, pp. 1039β1049, 2003. 19. L. D. Mosgaard, K. a. Zecchi, and T. Heimburg, Mechano-capacitive properties of polarized membranes,
Soft Matter , vol. 11, no. 40, pp. 7899β7910, 2015. 20. P.-C. Zhang, A. M. Keleshian, and F. Sachs, Voltage-induced membrane movement,
Nature , vol. 413, no. 6854, pp. 428β432, 2001. 21. I. Cinelli, M. Destrade, M. Duffy, and P. Mchugh, Electro-thermal equivalent 3D Finite Element Model of a Single Neuron,
IEEE Trans. Biomed. Eng. , vol. 9294, no. 99, 2017. 22. I. Cinelli, M. Destrade, M. Duffy, and P. Mchugh, Neurotrauma Evaluation in a 3D Electro-Mechanical Model of a Nerve Bundle, in , 2017, no. 1, pp. 2β5. 23. I. Cinelli, M. Destrade, M. Duffy, and P. Mchugh, Electro-Mechanical Response of a 3D Nerve Bundle Model to Mechanical Loads Leading to Axonal Injury, in , 2017, pp. 1β4. 24. A. L. Hodgkin and A. F. Huxley, A Quantitative Description Of Membrane Current And Its Application To Conduction And Excitation In Nerve,
J. Physiol. , vol. 117, pp. 500β544, 1952. 25. A. El Hady and B. B. Machta, Mechanical surface waves accompany action potential propagation,
Nat. Commun. , vol. 6, no. 6697, p. 6697, 2015. 26. I. Cinelli, M. Duffy, and P. E. McHugh, Thermo-electrical equivalents for simulating the electro-mechanical behavior of biological tissue, in , 2015, pp. 3969β3972. 27. C. Demerens et al. , Induction of myelination in the central nervous system by electrical activity., Proc. Natl. Acad. Sci. U. S. A. , vol. 93, no. 18, pp. 9887β9892, 1996. 28. P. D. Einziger, L. M. Livshitz, a. Dolgin, and J. Mizrahi, Novel cable equation model for myelinated nerve fiber,
First Int. IEEE EMBS Conf. Neural Eng. 2003. Conf. Proceedings. , vol. 52, no. 10, pp. 1632β1642, 2003. 29. S. A. Tawfik, D. Stefan Dancila, and E. Armanios, Unsymmetric composite laminates morphing via piezoelectric actuators,
Compos. Part A Appl. Sci. Manuf. , vol. 42, no. 7, pp. 748β756, 2011. 30. M. Staworko and T. Uhl, Modeling and simulation of piezoelectric elements - comparison of available methods and tools,
Mechanics , vol. 27, no. 4, pp. 161β171, 2008. 31. O. Alvarez and R. Latorre, Voltage-dependent capacitance in lipid bilayers made from monolayers.,
Biophys. J. , vol. 21, no. 1, pp. 1β17, 1978. 32. J. a Perge, J. E. Niven, E. Mugnaini, V. Balasubramanian, and P. Sterling, Why do axons differ in caliber?,
J. Neurosci. , vol. 32, no. 2, pp. 626β638, 2012. 33. I. Cinelli, M. Duffy, and P. E. McHugh, Thermo-electrical equivalents for simulating the electro-mechanical behavior of biological tissue, in
Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS , 2015, pp. 3969β3972. 34. B. Tahayori, H. Meffin, S. Dokos, A. N. Burkitt, and D. B. Grayden, Modeling extracellular electrical stimulation: II. Computational validation and numerical results.,
J. Neural Eng. , vol. 9, no. 6, p. 65006, Dec. 2012. 35. H. Meffin, B. Tahayori, D. B. Grayden, and A. N. Burkitt, Modeling extracellular electrical stimulation: I. Derivation and interpretation of neurite equations.,
J. Neural Eng. , vol. 9, no. 6, p. 65005, Dec. 2012. 36. J. Platkiewicz and R. Brette, A threshold equation for action potential initiation,
PLoS Comput. Biol. , vol. 6, no. 7, p. 25, 2010. 37. P. Boucher, B. JoΓ³s, and C. E. Morris, Coupled left-shift of Nav channelsβ: modeling the Na + -loading and dysfunctional excitability of damaged axons, pp. 301β319, 2012. 38. C. Cao, A. Yu, and Q.-H. Qin, Evaluation of effective thermal conductivity of fiber-reinforced composites,