Energy based mechano-electrophysiological model of CNS damage at the tissue scale

Injuries to the central nervous system (CNS) have become a major health challenge in our society. Such injuries can affect the electrophysiological properties, and thus functional behaviour of the nervous system. Most of the related computational studies to date have mainly been focussed on the cellular, subcellular or axonal scales. Although essential for understanding the characteristics of spinal cord injury and traumatic brain injury at the cellular level, these approaches suffer from a lack of scalability to meso- and macroscopic scenarios. To this end, this work presents a tissue level framework describing the coupled mechanical and electrophysiological behaviour of CNS tissue. The proposed mechanical large deformation constitutive model accounts for strain rate dependency, viscous effects and mechanical damage within a thermodynamically consistent framework. Coupled to it, a FitzHugh–Nagumo based model that includes axonal ﬁbre informed anisotropic conduc- tion describes the electrophysiological behaviour. The mechanical model is ﬁrst calibrated against experimental deformation measurements of spinal cord white matter. The com- plete framework is then calibrated and used to study the mechano-electrophysiological coupling in spinal cord samples subjected to tension-free relaxation tests by considering electrophysiological damage as a function of mechanical energetic terms. The ﬂexibility of the proposed model is illustrated with the case of mechano-electrophysiological cou- pling in spinal cord subjected to blast loading. After validation of the coupled model, the mechanical parameters are identiﬁed for white matter and an analysis of the inﬂuence of axonal dispersion on anisotropic AP conduction in an in silico rat model is presented. This work provides a novel mechano-electrophysiological framework able to model the mechanics of spinal cord and brain tissue and its translation into electrophysiological dysfunction through energetic terms at the tissue level. reserved.


a b s t r a c t
Injuries to the central nervous system (CNS) have become a major health challenge in our society. Such injuries can affect the electrophysiological properties, and thus functional behaviour of the nervous system. Most of the related computational studies to date have mainly been focussed on the cellular, subcellular or axonal scales. Although essential for understanding the characteristics of spinal cord injury and traumatic brain injury at the cellular level, these approaches suffer from a lack of scalability to meso-and macroscopic scenarios. To this end, this work presents a tissue level framework describing the coupled mechanical and electrophysiological behaviour of CNS tissue. The proposed mechanical large deformation constitutive model accounts for strain rate dependency, viscous effects and mechanical damage within a thermodynamically consistent framework. Coupled to it, a FitzHugh-Nagumo based model that includes axonal fibre informed anisotropic conduction describes the electrophysiological behaviour. The mechanical model is first calibrated against experimental deformation measurements of spinal cord white matter. The complete framework is then calibrated and used to study the mechano-electrophysiological coupling in spinal cord samples subjected to tension-free relaxation tests by considering electrophysiological damage as a function of mechanical energetic terms. The flexibility of the proposed model is illustrated with the case of mechano-electrophysiological coupling in spinal cord subjected to blast loading. After validation of the coupled model, the mechanical parameters are identified for white matter and an analysis of the influence of axonal dispersion on anisotropic AP conduction in an in silico rat model is presented. This work provides a novel mechano-electrophysiological framework able to model the mechanics of spinal cord and brain tissue and its translation into electrophysiological dysfunction through energetic terms at the tissue level.

Introduction
Injuries in the central nervous system (CNS) and, in particular, spinal cord injury (SCI) and traumatic brain injury (TBI) remain as an unsolved challenge with a significant socioeconomic impact on our society ( Sekhon and Fehlings, 2001;Jackson et al., 2004 ). SCI usually arises from mechanical loading leading to non-physiological axonal deformations within the spinal cord tissue. Axonal damage, especially at high strain rates, can in turn lead to electrophysiological alterations and, subsequently, partial or complete loss of sensory and motor function ( Anderson and Hall, 1993;Honmou and Young, 1995 ).
High deformation rates can be observed in many different loading scenarios such as during car accidents, falls and sports activities ( Maikos et al., 2008 ) or during blast loading scenarios (e.g., explosions) ( Warden et al., 2006 ). Therefore the characterisation of the mechano-electrophysiological coupling in CNS components is of great interest to evaluate the consequences of mechanical deformation in terms of functional and cognitive deficits.
To this end, experiments conducting longitudinal (uniaxial) strain on spinal cord specimens have been performed by Blight and DeCrescito (1986), Shi and Pyor (2002) and Shi and Whitebone (2006) . In all three effort s, axonal injury was evaluated by measurements of membrane permeability and reductions in the compound action potential (CAP), observing dependencies on both strain and strain rate. LaPlaca et al. (20 05,20 06) studied the influence of shear deformation at high rates in neural cells. These results showed that higher rates and shear strains result in higher neurites loss and membrane permeability. Morrison III and co-authors ( 20 06, 20 07, 20 09a, 20 09b, 2015 ) carried out biaxial deformation experiments on brain slices of different regions. A direct relationship between strain and cell death was found for all the brain regions tested and strain rate was identified to play a significant role in some regions such as cortex but not in others such as hippocampus ( Elkin and Morrison III, 2007 ). In addition, this group developed a methodology to evaluate alterations in hippocampal network activity due to mechanical deformations ( Kang et al., 2015 ). Electrophysiological alterations have also been studied within the context of blast loading ( Connell et al., 2011a, b ). These experiments showed that blast exposure results in deformations at very high rates leading to electrophysiological alterations. More recently, Tagge et al. (2018) conducted unilateral impacts on mice and observed bilateral impairments in axonal conduction velocity and electrophysiological deficits in hippocampus and prefrontal cortex.
Numerical models have also been explored to complement such studies with data not necessarily directly available through experimental means. In particular, SCI finite element (FE) models have been developed for the translation of knowledge from animal SCI studies to humans ( Jannesar et al., 2018 ). In addition, they provide essential tools to shed light in the quantitative evaluation of injury to the spinal cord ( Karimi et al., 2017 ). These models are usually focussed on mechanical aspects and their consequences on the CNS ( Kleiven and Hardy, 2002;Wright et al., 2013;Giordano et al., 2017;Alvarez and Kleiven, 2018 ) and, eventually, on functional and cognitive deficits . However, there are only few effort s that model the direct influence of mechanics on the electrophysiology, and most of them are mainly focussed on the cellular, subcellular or axonal scales. One of the first attempts to model this coupled behaviour was done by Boucher et al. (2012) , who analysed the left shift dynamics of Na v channels induced by trauma, and by Babbs and Shi (2013) , who simulated mild retraction of myelin due to stretch and crush injuries. Nevertheless, these authors did not directly couple the electrophysiological alterations to mechanical variables but rather modified the electrophysiological parameters indirectly as a function of the new mechanical state. To this end, Jerusalem and co-authors ( 2014 and 2015 ) developed a 1D finite difference framework to model the effects of strain and strain rate on axons under uniaxial loading. This model was extended to deal with 3D non-axisymmetric mechanical deformations within a FE framework ( Kwong et al., 2018 ). Similar effort s developing 3D non-axisymmetric coupled models have been recently proposed ( Cinelli et al., 2018a, b, c ). Finally, Ng et al. (2017) linked the mechanics of FE models at the macroscopic level to neurologic injury at the cell level while  linked it directly to medium-term cognitive deficits through geometrical correlation of oxidative stress regions.
Although essential for understanding SCI and TBI at the cellular level (or directly at the cognitive scale in the latter reference), these approaches suffer from a lack of applicability to tissue scale electrophysiology scenarios. To this end, we present a tissue level framework to describe the coupled mechanical and electrophysiological behaviour of axonal based tissues. This framework combines a large deformation mechanical constitutive model with an electrophysiological tissue model. The mechanical part accounts for strain rate dependency, viscous effects and continuum mechanical damage and is developed within a thermodynamically consistent framework. The electrophysiological behaviour is described by a FitzHugh-Nagumo (FHN) based model ( FitzHugh, 1961;Nagumo et al., 1962 ) that includes anisotropic conduction following fibre orientation. The mechanical model is first calibrated against experimental measurements of the stress-strain response of spinal cord tissue under a wide range of strain rate conditions, and the electrophysiological parameters are obtained from the literature. Then, the complete framework is calibrated (for the coupling parameters) and used to analyse the mechanoelectrophysiological coupling in spinal cord samples subjected to tension-free relaxation tests by considering electrophysiological damage as a function of mechanical energetic terms. The predictive capability of the framework is illustrated with the complex scenario of mechano-electrophysiological coupling in blast-loaded spinal cord. Finally, the mechanical parameters are identified for white matter and an analysis of the influence of axonal dispersion on anisotropic AP conduction is presented.

Mechano-electrophysiological constitutive framework
A 3D mechano-electrophysiological constitutive model is proposed for finite deformations and anisotropic diffusion. The mechanical formulation is developed within a thermodynamically consistent framework whose constitutive relations derive from energy potentials. The electrophysiological formulation is based on the FHN model ( FitzHugh, 1961;Nagumo et al., 1962 ) and incorporates anisotropic electrophysiology defined by the axonal-based preferred orientations. Finally, the coupling between mechanics and electrophysiology is introduced through energetic terms.

Mechanical constitutive framework
In this section, we describe the mechanical framework and the corresponding constitutive relations to define the mechanical behaviour of CNS tissue. To this end and in order to facilitate the implementation of the quasi-incompressibility condition present in such tissues, the deformation gradient is split into a volume-changing component, F vol = J 1 / 3 I , and an isochoric component, F iso = J −1 / 3 F : where J = det F is the Jacobian and I is the second-order identity tensor. Therefore, a dilated intermediate configuration ¯ can be reached from the reference configuration o by uniquely accounting for the mechanically induced dilation of the continuum body, see Fig. 1 a. Moreover, the overall response of the material is understood as the summation of a purely hyperelastic contribution, denoted as the Equilibrium (Eq) branch, and a visco-hyperelastic contribution, denoted as the Non-Equilibrium (NEq) branch. The viscous resistance is defined by a hyperelastic spring in series with a nonlinear dashpot, see Fig. 1 b. According to this rheological representation of the framework, the isochoric deformation gradient is further decomposed into elastic ( F e ) and viscous ( F v ) components as By virtue of the multiplicative decomposition of F iso , a second intermediate configuration is defined. This configuration is referred to as dilated relaxed configuration ¯ and accounts for both mechanical-induced dilation and viscous deformation. Therefore, the spatial velocity gradient L can be written as where L vol and L iso are its volumetric and isochoric components, respectively. The isochoric component can be further decomposed as where L e is the elastic velocity gradient and . In this work, ¯ is assumed invariant to rigid body rotations of the current configuration , that is Boyce et al., 1988 ).

Thermodynamics
The mechanical framework is developed to be thermodynamically consistent with the aim of coupling the mechanical behaviour of the soft tissue with its electrophysiological behaviour through energetic terms. To this end, we postulate the existence of a Helmholtz free energy per unit reference volume = e − T η, where e is the specific internal energy in o , T is the absolute temperature and η is the specific entropy in o . This energy function is decomposed into a volumetric component vol , an isochoric component related to the Equilibrium component Eq , an isochoric component related to the Non-Equilibrium component NEq and a thermal component t as The isochoric components of the energy are weighted by the reduction factor ( 1 − D ) ( Kachanov, 1958;Simo, 1987 ), were D ∈ [0, 1] is an internal variable associated with a continuous material damage or degradation due to mechanical deformation. Note that the Non-Equilibrium component can be alternatively defined in terms of Although the thermal component of the free energy is neglected in this work, it is added not only for the sake of thermodynamics consistency but also to allow for the definition of the specific entropy. In this regard, the proposed framework allows for the extension of the model to recent theories that, instead of considering the classical view of ionic currents determining the action potential (AP) propagation, understand the electrophysiological signal as an adiabatic pulse whose physical interpretation is strongly related to the entropy evolution ( Heimburg and Jackson, 2005 ). Moreover, the decomposition of the free energy function into volumetric and isochoric components is chosen to facilitate the choice of volumetric or isochoric components as governing variables of the electrophysiological alteration in an independent manner.
Taking the Helmholtz free energy definition from Eq. (5) along with the first and second Principles of Thermodynamics, the following expression for the Clausius-Plank inequality is reached where P is the first Piola-Kirchhoff stress and P = I − 1 3 F −T F defines a fourth-order projection tensor with F operating as a metric tensor, and I denotes a fourth-order unity tensor. The inequality presented above must holds for any arbitrary deformation. In this regard, the first and second terms satisfy this requirement by applying the standard Coleman-Noll procedure ( Coleman and Noll, 1963;Coleman and Gurtin, 1967 ) to obtain the following constitutive relations for the stress state and the specific entropy per unit volume: The first Piola-Kirchhoff can be further decomposed into its volumetric and isochoric terms as ∂J denotes the hydrostatic pressure. Moreover, the third and fourth terms of Eq. (6) refer to the inelastic dissipation due to viscous effects associated with the Non-Equilibrium thermodynamics ( D v int ) and to the irreversible damage dissipation ( D d int ), respectively. Eqs.
(1 -12 ) define the general mechanical framework. To complete the constitutive description of this model, the specification of the energy potentials chosen as well as the definition of the viscous and damage flow rules are needed.

Helmholtz free energy functions
Here, specific forms of the Helmholtz free energy functions vol , Eq , NEq are provided. The volumetric part of the strain energy function is defined by use of the bulk modulus of the material K : Regarding the energy terms associated to the isochoric components, the Ogden energy function ( Ogden, 1972 ) is used as it has been shown to accurately describe the mechanical behaviour of CNS soft tissues ( Garcia-Gonzalez et al., 2018 ) where λ iso i and λ e i are, respectively, the principal stretches of the isochoric and elastic parts of the deformation gradient with i ∈ {1, 2, 3}. In the Ogden model, the equivalent shear moduli μ i α i /2 of each mechanical branch is given by the material parameters μ i and α i , with i ∈ { Eq, NEq }. Note that the physical consistent condition μ i α i ≥ 0 must be verified ( Ogden et al., 2004 ). Also, since the formulation is written in terms of principal stretches, the expressions for the First Piola-Kirchhoff stress tensors ( Eqs. (10.2) and ( 10.3 )) can be given as where N iso i and N e i are the corresponding eigenvectors of the isochoric and elastic components expressed in the reference configuration, respectively.
The thermal energy term t ( T ) is associated with the entropy evolution and is neither defined nor used in this work.

Viscous. flow rule
The viscous stretch rate D v is defined as where ˙ γ v is the visco-elastic multiplier and N v is a tensor that provides the direction of the viscoelastic flow in the deformed configuration : where σ dev NEq = ( P NEq F T ) dev is the deviatoric Cauchy stress tensor associated with the non-equilibrium thermodynamics and σ dev NEq = tr ( σ dev NEq σ dev NEq ) . The visco-elastic multiplier ˙ γ v describes the effective flow rate and is constitutively defined as ( Prevost et al., 2011 ) where G base is a dimensional scaling constant, 0 < α 1 is a constant incorporated to eliminate singularity and n and c are material constants.

Damage evolution
The flow rule of the damage variable D must be provided to complete the mechanical framework. To this end, we follow Simo (1987) and define an internal variable governing the damage evolution as This variable is used to define the damage criterion ϕ that defines the damage surface condition: where max is the maximum value between a critical energy threshold crit and the maximum ( t ) reached over the past deformation history up to the current time t : Finally, the evolution of the damage variable D is defined as Note that this formulation ensures that damage grows linearly as 1/ s as soon as reaches crit , where s is a saturation energy threshold, and full damage is asymptotically reached as the internal energy goes towards infinity as Although the specific energies per unit volume crit and s could loosely represent energetic thresholds for the onset of permanent damage of the micromechanical components of the CNS tissue, and its complete damage, respectively, it is worth emphasising that the evolution of the mechanical damage as defined here is phenomenological in nature to bypass purposefully the need for micromechanical consideration.

Electrophysiological constitutive framework
The electrophysiological behaviour of axons and nerve fibres is usually defined by different variations of the Hodgkin-Huxley (HH) model ( Jerusalem et al., 2014;Garcia-Grajales et al., 2015 ). Such models take into account features of the different ionic currents and gate variables to explain the propagation of the AP. However, especially in cardiac mechanics, the modelling of the electrophysiological behaviour at tissue level has often been addressed by simplifications of the HH model that reduce the phenomenological variables associated to different ion currents to one. In this regard, the most popular approaches to model AP propagation at tissue level are the Aliev-Panfilov model ( Aliev and Panfilov, 1996;Costabal et al., 2017 ) and FHN model ( Engelbrecht et al., 2018 ). In this work, we use a 3D diffusion model based on FHN as where φ is the normalised voltage, ∇ X is the material gradient with respect to the reference configuration o , D is the conductivity tensor and R φ is the material ionic current, both expressed in in the reference configuration o . The term R φ is defined according to FitzHugh (1961) and Heimburg and Jackson (2005) as where 0 < a < 1 is an activation parameter, I stim is a stimulus current and r is the recovery current following: where ɛ is the time-scale difference and 0 < b < 1 is an activation parameter. Finally, in order to include anisotropic AP conduction, we define the conductivity tensor as where D iso and D ani are, respectively, the isotropic and anisotropic contributions of the material conductivity; and A axon is the axonal structure tensor. The axonal structure tensor is defined following the formulation used by Garcia-Gonzalez et al. : Here, ξ is a dispersion parameter that depends on the fractional anisotropy FA and a axon is the preferred axonal orientation. When ξ adopts a value of 1/3, an isotropic distribution of the axons is expected and, in this case, the structure tensor is spherical. When ξ adopts a value of 0, an ideal coalignment of the axons is expected and the structure tensor reduces to A axon = a axon a axon . Note that this structure tensor approach can be also used to include anisotropy in the mechanical response due to axonal orientation and dispersion; for more details, see  .

Mechano-electrophysiological coupling
To complete the mechano-electrophysiological framework, energy-based terms extracted from the mechanical framework are incorporated into the electrophysiological formulation to induce a continuous reduction in the AP due to mechanical deformation. To this end, we postulate that the reduction in AP propagation is related to the mechanical damage D and an energy-based variable, hereafter referred to as compound energy ω, defined by a phenomenological function dependent on the total isochoric free energy ( Eq + NEq ) and the inelastic dissipation term D v int derived in the Clausius-Plank inequality, see Eq. (6) . Therefore, the compound energy variable can be expressed as ω = ˆ ω ( Eq , NEq , D v int ) and Eq. (25) can be rewritten as where 0 < ˆ b < 1 is the activation parameter depending on the mechanical damage and energetic terms from the mechanics. Based on the power-like behaviour observed in the subsequent combined experimental-numerical analysis, this functional dependence was chosen to adopt the following form where b o is the activation parameter in the healthy reference and p 1 , p 2 and p 3 are material parameters. Note that when the mechanical damage parameter D tends to 1 (tissue fully damaged), the activation parameter ˆ b tends to infinity, impeding the propagation of the electrophysiological signal.
To complete the framework, the function ω = ˆ ω ( Eq , NEq , D v int ) must be defined from the combination of experimental data and numerical simulations (analysed in the following section).

Mechanical calibration
The calibration of the mechanical response of spinal cord tissue is based on characterisation experiments published by Fradet et al. (2016) . These authors conducted uniaxial compression tests on spinal cord specimens under different loading conditions covering strain rates from 0.5 s −1 to 50 s −1 , see Fig. 2 . The experimental results show large deformation of the tissue and nonlinear mechanical behaviour. In addition, there is a clear strain rate dependency: higher strain rates result in tissue stiffening. Moreover, at high deformations (strains larger than 55%) the tissue experiences a continuous degradation in its mechanical response, which also depends on strain rate and is associated to mechanical damage. To capture the large deformation and nonlinear response, the Ogden free energy function defined in Eqs. (14.1) and ( 14.2 ) were found to faithfully describe the stress-strain shape of the curves, and thus adopted here. The strain rate dependency is captured by the nonequilibrium branch through Eqs. (16) -(18) , and the continuous material degradation is captured by the continuum damage model through Eqs. (19) -(22) . In addition, the continuum damage model predicts adequately the strain rate dependency on damage onset through the consideration of the energy term NEq ( F e ) in the driving force. It is worth mentioning that the model does not capture the initial s-shape region of the stress-strain curve. However, as the mechano-electrophysiological coupling is defined in terms of free energy, viscous dissipation and mechanical damage, the coupling is the weakest in this region where, therefore, the stress-strain predicting accuracy is less relevant for the electrophysiological prediction.
The identification of the mechanical model parameters is achieved by making use of a sequential quadratic programming optimisation method ( fmincon , in MATLAB terminology). The resulting mechanical parameters are provided in Table 1 , and the model predictions are shown together with the experimental curves in Fig. 2 . Note that the material parameters for the Ogden free energy ( μ i and α i ) are both negative and satisfy the physical consistent condition μ i α i ≥ 0, and are consistent

Table 2
Material parameters for the electrophysiological behavior of spinal cord tissue identified from the work published by Engelbrecht et al. (2018) . with experimental observations on axonal based tissues that show stiffer material response in compression with respect to tension ( Budday et al., 2017a,b ).

Electrophysiological calibration
The electrophysiological calibration requires the identification of the FHN model parameters a , ɛ and b o ; as well as the conductivity parameters. With the aim of evaluating the consistency of the proposed approach and identifying the model parameters, the electrophysiological part of the 1D model proposed by Engelbrecht et al. (2018) for nerve fibres is first implemented in a finite difference framework. An electrophysiological simulation of the AP propagation is then conducted and compared with the results of the FE framework proposed in this work when assuming the axonal orientation along the longitudinal direction of the spinal cord tissue ( a axon ) and a FA value of one (fibres perfectly aligned). In addition, the equivalent conductivity used by Engelbrecht et al. (2018) is assumed to be the sum of both isotropic and anisotropic conductivity components as D total = D iso + D ani . Fig. 3 shows the matching prediction of the electrophysiological response of nerve fibres by the model proposed in the previous section ( Eqs. (23 -28 )) against the model proposed by Engelbrecht et al. (2018) for such tissues. The model parameters used are provided in Table 2 . Note that all the equations are expressed in dimensionless terms and a mapping of the resulting fields and source terms can be straightforwardly obtained to translate these results into a realistic physiological domain, see Costabal et al. (2017) for details. Time and spatial convergences were verified (see Section 4.1.1 ).

Mechano-electrophysiological coupling calibration
This section evaluates the predictive capability of the presented mechano-electrophysiological framework and analyses the influence of mechanical deformation on the electrophysiological response in terms of AP propagation. To this end, the experimental work done by Shi and Whitebone (2006) is taken as the basis for this study. In this work, the authors developed an experimental model to study the functional and structural significance of tensile strain injury on guinea pig spinal cord white matter. Shi and Whitebone (2006) provided the CAP for each healthy spinal cord specimen and its reduction immediately after the tensile test and continuously for the subsequent 30 min. Six different loading conditions are studied combining strain magnitude and rate: strain magnitudes of 0.25 (mild), 0.50 (moderate) and 1.00 (severe); slow strain rates (0.0 06-0.0 08 s −1 ) and high strain rates (355-519 s −1 ). Table 3 Material parameters for the mechano-electrophysiological coupling of guinea pig spinal cord tissue as calibrated against the experiments reported by Shi and Whitebone (2006) .

FHN Activation Parameter
The experiments conducted by Shi and Whitebone (2006) show that during loading conditions, the CAP experiences a reduction in amplitude as the strain, strain rate or both increase. The loading step is followed by a recovery of CAP amplitude during the unloading step, eventually leading to a permanent amplitude reduction after complete relaxation. To explain these tendencies, we postulate that the mechano-electrophysiological coupling depends on a recoverable component and a permanent (or unrecoverable) component. The former is assumed to depend on the total isochoric free energy ( Eq + NEq ) while the latter is assumed to depend on the inelastic dissipation ( D v int ) and the mechanical damage of the tissue ( D ). The dependency of the recoverable component on the isochoric energy implies a continuous reduction of the electrophysiological response under loading conditions and a continuous recovery during unloading conditions. In addition, it accounts for strain rate dependency through the contribution of the non-equilibrium branch, providing rate-dependent predictions of the CAP reduction in agreement with experiments: larger CAP reductions for higher strain rates. Regarding the unrecoverable component, we postulate that only the inelastic dissipation that takes place during loading conditions affects the electrophysiological response of the tissue. This means that dissipative energy terms arising from non-loading tissue relaxation (negative isochoric energy rates) do not contribute to variations in electrophysiology. This assumption is introduced to reproduce the experimental observations that further tissue damage does not a priori occur during free unloading. Finally, the unrecoverable component also depends on the level of material degradation through the variable D , representing the CAP reduction associated with tissue damage resulting in an increase of permeability ( Shi and Whitebone, 2006;Laplaca et al., 2018 ). These assumptions are taken together to define the compound energy ω as where H( ˙ Eq + ˙ NEq ) is the Heaviside unit step function related to the isochoric energy rate. Note that with the introduction of this term, the inelastic dissipation energy is accounted for only under loading conditions and not under free unloading or relaxation of the tissue. One has: Combining Eqs. (30 -33 ), the mechano-electrophysiological framework is completely defined, pending the calibration of the coupling parameters p 1 , p 2 and p 3 so as to reproduce the multiphysics experiments carried out by Shi and Whitebone (2006) . The final values for the coupling parameters are obtained through the same optimisation process as in the mechanical calibration (sequential quadratic programming optimisation method) and are provided in Table 3 . The model predictions for the CAP reduction, depending on the compound energy computed for the six different loading conditions and under stressed and relaxed states, are plotted together with the experimental results in Fig. 4 . Note that, following the methodology proposed to identify the different parameters, the model can easily be calibrated for other tissues such as white matter or grey matter. In addition, note that the CAP represents the sum of individual AP within the spinal cord tissue. Therefore, the alteration of the electrophysiological response must be understood as a structural damage wheri, some of the axons within the tissue loose their availability to fire an AP, resulting in a reduction in the overall propagation of the electrical pulse.

Validation: mechano-electrophysiological response of spinal cord to blast loading
The mechano-electrophysiological framework developed in Section 2 and calibrated in Section 3 was used to analyse the coupling between mechanical deformation and electrophysiological damage under uniaxial tensile loading. Here, this framework is used to analyse more complex scenarios and to validate its prediction capability for high rate loading. More specifically, we study the effects of blast loading on spinal cord tissue and the subsequent reduction in CAP. To this end, the experimental work by Connell et al. (2011a) is taken as the basis for this study. In this work, the authors developed an experimental model to study the functional and structural significance of blast-induced deformation on guinea pig spinal cord white matter. Connell et al. (2011a) exposed isolated spinal cord specimens to three blast wave intensities: 23 kPa (mild), 41 kPa (moderate) and 65 kPa (severe); and provided the CAP for each healthy spinal cord and its reduction during the subsequent post-loading 30 min.

Numerical model
The blast experiments are simulated with a combined Lagrangian-Eulerian 3D FE model developed in Abaqus/Explicit ( Abaqus v6.14, 2014 ). The model consists of a Lagrangian mesh that represents the spinal cord specimen, a rigid body that represents the base where the specimen lays on and a Eulerian mesh representing the mass of air through which the shock wave propagates and impacts the spinal cord, see Fig. 5 . The geometrical features of this set up are chosen accordingly to the details given by Connell et al. (2011a) .
Here, the quasi-incompressible behaviour of spinal cord is introduced in the FE model by assuming a Poisson's ratio of ν = 0 . 49 ( Maikos et al., 2008 ). The equivalent shear modulus at high deformation rates can be approximated as μ = μ Eq α Eq / 2 + μ NEq α NEq / 2 and the bulk modulus used for this simulation can be calculated as: The remaining parameters used in these simulations for both mechanics and electrophysiology are identified in the previous section. A convergence analysis was conducted to evaluate the time and spatial discretisation. In this regard, a mesh size of 0.2 mm was found to ensure spatial convergence for both mechanical and electrophysiological convergence (see Section 3.2 ). The critical time step used for the implicit integration of the electrophysiological simulation after the blast is 0.005 ms.

Results
In order to validate the mechanics of the spinal cord tissue when subjected to blast, the experimental measurement of an overall deformation of the tissue versus time is compared with the deformation predicted by the proposed model, see Fig. 6 for the severe loading case.
For this case, the maximum rate occurs within the first 80-120 μs for both the experiments and numerical predictions, and the maximum deformation magnitude reached is around 60% with a predominant compression component. The measurement of the deformation is computed following the experimental methodology used in the experiments, see Connell at al. (2011a) . The deformation process is spread over a longer period in the numerical predictions than in the experimental measurements. However, the strain rate remains of the same order of magnitude for both cases. The mechanical deformation of the tissue under the three different blast exposure levels after 30 min results in a distribution of compound energies within the impacted tissue region with a mean value of ∼50 J/m 3 for mild conditions, ∼80 J/m 3 for moderate conditions and ∼300 J/m 3 for severe conditions. Finally, when these loading conditions are simulated with the complete mechanoelectrophysiological model, the mechanical deformations from the blast exposures result in long term (30 min) CAP reductions of 21.84% for mild conditions, 35.79% for moderate conditions and 54.82% for severe conditions. These results are compared with the corresponding experimental values reported by Connell et al. (2011a) in Fig. 7 , showing a good agreement between them and verifying the subsequent good predictive capability of the framework proposed.

Analysis of anisotropic AP conduction: influence of axonal dispersion
This section presents a qualitative analysis of the influence of axonal dispersion on the anisotropic AP conduction. These simulations must be understood as a first approach to model the electrophysiological behaviour in CNS tissues. For a more accurate modelling of this response, a connectivity tensor is needed on top of the proposed formulation to determine discrete paths for the AP propagation. To this end, a rat hippocampus brain slice, similar to the organotypic hippocampal slice culture studied experimentally by Morrison et al. (2006 ) is chosen here. The in silico model of the slice is developed from magnetic resonance imaging and diffusion magnetic resonance imaging (dMRI) . It is worth emphasising that the slice chosen in this work actually corresponds to rat grey matter. However, as diffusion tensor imaging data was captured for this region by  , white matter material properties were used as a first approximation. Doing so, a mean axonal orientation a axon is defined for each element to compute the local structure tensors A axon . The fractional anisotropy FA is then varied to study the influence of axonal dispersion on the AP conduction.
The electrophysiological parameters used in these simulations are the same set employed in the previous sections (see Table 2 ) and D iso is taken to be zero in Eq. (26)  mechano-electrophysiological experiments that could allow for the calibration of the coupling parameters, the same set of parameters used for the white matter spinal cord is used here (see Table 3 ). The mechanical parameters are calibrated against quasi-static and very high rate compression tests conducted by Pervin and Chen (2009) on white matter samples. The model predicts the nonlinear mechanical behaviour of the tissue and the strong strain rate sensitivity, see Fig. 8 for comparison with experimental data. However, as there is no evidence of mechanical damage in these experiments, those were not accounted for in this simulation. The calibrated parameters can be found in Table 4 . The brain slice of the rat is then subjected to biaxial loading representative of the deformation process during impacts ( Morrison et al., 2006 ). In this regard, an equivalent strain of 21% is imposed biaxially to the brain tissue at a strain rate of 10 0 0 s −1 . Note that these values are in the expected range for TBI injury ( Wright et al., 2013 ). Finally, the resulting model is used to simulate the AP conduction before and after the loading process by triggering an electrical signal in an arbitrary location of the slice. The simulations are run for perfectly alignment of axons within the elements (FA = 1), and for fully disperse axons within the elements (FA = 0). Fig. 9 shows the AP conduction for both undeformed and deformed conditions, and the influence of the axonal dispersion on this conduction. High FA values result in a strong anisotropic AP conduction while low FA values result in an isotropic AP conduction. The deformed results suggest that regions with dispersed axons  are more likely to suffer an alteration in the electrophysiological response than regions where the axons are well aligned. This observation can be explained by the fact that diffuse propagation of the electrical signal for low FA values offer more signal "rerouting" opportunities, and hence signal alterations, than its high FA counterpart. As a conclusion, this illustrative example demonstrates that higher energies (i.e., higher deformation and/or deformation rates) are needed to significantly alter the electrophysiological response in regions where the axons are well aligned.

Discussion
The proposed model is calibrated against tension-relaxation tests at different deformation levels and rates under certain assumptions that couple mechanical deformation with electrophysiological alterations. The mechanically-induced electrophysiological deficit is understood as the combination of a recoverable component, defined here as the total isochoric free energy, and a permanent component, defined here as dependent on viscous dissipation and mechanical damage. These variables are derived from energetic terms per unit volume, allowing to homogenise traditional approaches from the cell scale to the tissue scale. Furthermore, the proposed model takes into account strain rate sensitivity in the evolution of the different energetic variables (and mechanical damage) allowing for the prediction of electrophysiological alterations within a wide range of both quasi-static and dynamic deformation processes. The flexibility of the model is demonstrated against complex mechanical loading, i.e., blast loading, for which the deformation rate is very high, where the coupled model is able to predict the altered electrophysiological response in good agreement with experiments. Moreover, the model takes into account the influence of axonal orientation and dispersion on the propagation of the AP within the tissue.
Thanks to its continuum nature, the proposed mechano-electrophysiological framework is particularly relevant to clinical problems that must be addressed through a tissue-scale approach, i.e., bypassing all microstructural considerations. Among the clinical problems that involve low deformation rates, the proposed model can be used, for example, together with experiments to elucidate the causes of carpal tunnel syndrome. Carpal tunnel syndrome is the result of the compression of the median nerve within the carpal tunnel of the wrist, that results in several symptoms such as pain, numbness and tingling in the fingers ( Coppieters and Alshami, 2007 ). Indeed, while a full consideration of all the axons involved in this disease is impossible, the proposed approach could allow for the simultaneous consideration of the wrist mechanical simulation along with the electrophysiology of the compressed nerve. Other dynamic scenarios that result in mechanical deformation of CNS tissue and require a meso-to macroscale approach include car or bike accidents, fall of a person or even blast injury. In this regard, a model was recently presented by the authors, relating the evolution of mechanical variables during dynamic deformation processes to cognitive deficits . The present mechano-electrophysiological model, as formulated at the tissue scale, could be implemented for the white matter incorporating the axonal information from dMRI to analyse potential alterations in the brain electrophysiological behaviour, eventually rationalising the observed cognitive deficits.

Conclusions
The main contribution of this work is the development of a mechano-electrophysiological framework to model the effects of mechanical deformation and tissue damage on the resulting electrophysiological deficits in central nervous system tissues. The mechanical formulation is built on the basis of a thermodynamically consistent framework for finite deformations that takes into account strain rate dependency, viscous effects and continuum mechanical damage. The electrophysiological formulation is defined by a FitzHugh-Nagumo based model that includes anisotropic conduction following fibre orientation. Energetic terms and mechanical damage, both derived from the mechanics, are included in the electrophysiological formulation to couple both physics and complete the framework.
Both mechanical and electrophysiological models are calibrated with experiments and published work available in the literature. The coupling between mechanics and electrophysiology is analysed by numerically reproducing published experiments of spinal cord tensile tests. The complete framework is validated against a complex blast-loading scenario. A good agreement between experimental results and model predictions is obtained for both tension-free relaxation and blast tests, showing a reduction in the CAP magnitude and propagation with increased strain and strain rate levels of deformation. Finally, the mechanical model is calibrated for white matter against experimental data and an analysis of the influence of axonal dispersion on anisotropic AP conduction in an in silico rat model is presented.
This novel mechano-electrophysiological framework eventually allows for the modelling of the mechanics of central nervous system tissues and its translation into electrophysiological dysfunction through energetic terms.