Computational model of the mechanoelectrophysiological coupling in axons with application to neuromodulation

For more than half a century, the action potential (AP) has been considered a purely electrical phenomenon. However, experimental observations of membrane deformations occurring during APs have revealed that this process also involves mechanical features. This discovery has recently fuelled a controversy on the real nature of APs: whether they are mechanical or electrical. In order to examine some of the modern hypotheses regarding APs, we propose here a coupled mechanoelectrophysiological membrane ﬁnite-element model for neuronal axons. The axon is modeled as an axisymmetric thin-wall cylindrical tube. The electrophysiology of the membrane is modeled using the classic Hodgkin-Huxley (H-H) equations for the Nodes of Ranvier or unmyelinated axons and the cable theory for the internodal regions, whereas the axonal mechanics is modeled by means of viscoelasticity theory. Membrane potential changes induce a strain gradient ﬁeld via reverse ﬂexoelectricity, whereas mechanical pulses result in an electrical self-polarization ﬁeld following the direct ﬂexoelectric effect, in turn inﬂuencing the membrane potential. Moreover, membrane deformation also alters the values of membrane capacitance and resistance in the H-H equation. These three effects serve as the fundamental coupling mechanisms between the APs and mechanical pulses in the model. A series of numerical studies was systematically conducted to investigate the consequences of interaction between the APs and mechanical waves on both myelinated and unmyelinated axons. Simulation results illustrate that the AP is always accompanied by an in-phase propagating membrane displacement of ≈ 1 nm, whereas mechanical pulses with enough magnitude can also trigger APs. The model demonstrates that mechanical vibrations, such as the ones arising from ultrasound stimulations, can either annihilate or enhance axonal electrophysiology depending on their respective directionality and frequency. It also shows that frequency of pulse repetition can also enhance signal propagation independently of the amplitude of the signal. This result not only reconciles the mechanical and electrical natures of the APs but also provides an explanation for the experimentally observed mechanoelectrophysiological phenomena in axons, especially in the context of ultrasound neuromodulation. DOI: 10.1103/PhysRevE.99.032406


I. INTRODUCTION
In the mid-20th century, Alan Hodgkin and Andrew Huxley published their distinguished Hodgkin-Huxley (H-H) model [1], eventually leading to the Nobel Prize in Physiology or Medicine a decade later. Fit against experimental data obtained from squid giant axons, the H-H model uses a set of nonlinear differential equations to capture the ionic mechanisms underlying the initiation and propagation of a nerve electrical pulse or action potential (AP). The H-H model has since become the research gold standard worldwide for the modeling of the electrical characteristics of neurons.
Since the late 1970s, however, experimental measurements of other various nonelectrical aspects occurring during APs have revealed that a nerve pulse involves more than just its electrical component [2][3][4][5][6][7]. In particular, the AP has been observed to be accompanied by a small and rapid surface * antoine.jerusalem@eng.ox.ac.uk motion [8][9][10][11][12][13][14][15][16][17][18][19][20] whose peak is in phase with that of the AP, with a magnitude in the nanometre range [21]. Marked longitudinal shrinkage of the nerve fiber has also been observed simultaneously with radial swelling [7]. Building on these experimental observations, a new perspective began to arise, essentially aimed at studying the phenomenon of AP propagation in nerve fibers from a thermodynamics standpoint [7,22]. In this approach, the phase transition of the lipid bilayers in nerve fibers plays a significant role during nerve pulse propagation [23]. Recently, the soliton hypothesis proposed by Thomas Heimburg and Andrew Jackson further reinforced the thermodynamics perspective [24]. The soliton hypothesis considers the AP as a localized adiabatic "mechanical density pulse" that maintains its shape during propagation with no attenuation. The key concept of the soliton hypothesis is the balancing of the pulse's dispersive property (frequencydependent wave speed) by its nonlinearity (due to the increase in compressibility as a result of membrane lateral compression) near the phase transition regime of the lipid bilayers [21,25]. Experimental support of the soliton model includes the measured temperature pulse and the absence of heat release during an AP. While the soliton model appears to disagree with the H-H model on the fundamental nature of the AP, the fact that the H-H model was essentially conceived as a phenomenological model for the electrical signal alone, along with the increasing body of evidence pointing toward the existence of mechanical pulses (MPs) in cell membrane, highlights the need to reconsider the AP as a full multiphysics problem [26][27][28].
Biological flexoelectricity, which includes both direct and reverse flexoelectricity, has been suggested as a key player in the coupling between mechanics and electrophysiology in biological cells [29]. Direct flexoelectricity is a unique property usually observed in liquid crystals and refers to the existence of a spontaneous electric polarization in dielectric materials submitted to an applied mechanical strain gradient [30]. Conversely, the reverse flexoelectric effect refers to a strain gradient (usually in the form of bending or distortion) induced by a substantial external electric field applied across the material [31,32]. In lipid bilayers and cell membranes, the phospholipid molecules which arrange themselves into two layers are charged with dipole moments. These fluidlike membranes can undergo a variety of mechanical deformations such as bending and thickness change due to external mechanical stimulation, resulting in a redistribution of the dipole moments and the development of a surface polarization [33]: the direct flexoelectric effect [34,35]. Reversibly, a substantially large external electric field across the membrane may also redistribute and reorientate these molecules due to electrical attraction and repulsion, resulting in a change in the membrane surface curvature: the reverse flexoelectric effect. Flexoelectricity is found to play a significant role in a number of biological functions, such as ion channels, thermal fluctuation, vesicle equilibrium, and hearing functions [10,29,[36][37][38][39][40]. In neurons in particular, it has been suggested that both the above-mentioned membrane movement during AP propagation, as well as the membrane polarization during an AP, could be related to the flexoelectricity [36,41].
Most recently, different numerical models have been proposed to study the coupling between the electrophysiological and mechanical properties of the membrane during AP initiation and propagation in both healthy and damaged neurons [27,28,[42][43][44][45][46][47][48]. However, to the best of the knowledge of the authors, none provided a numerical model able to investigate thoroughly the interactions between mechanical and electrical waves in the context of neuromodulation. To this end, we propose here a coupled mechanoelectrophysiological membrane finite-element model for neuronal axons. The model consists of four components: (1) an electrophysiological pulse model simulating the electrical AP propagation in axons, (2) a MP model simulating the propagation of deformation of axon membrane, (3) a reverse flexoelectric model linking membrane bending to the alteration of membrane potential, and (4) a direct flexoelectric model linking membrane polarization to the membrane mechanical strain gradient. Details of the modeling methodology, simulation results of pulse conduction and interaction on both unmyelinated and myelinated axons, as well as related discussions, are introduced in the following.

A. Overview
The geometry of the axon is modeled as a one-dimensional (1D) axisymmetric cylinder of initial length L, radius R, with a thin wall representing the axon membrane with a constant thickness H. A local cylindrical coordinate system (z, r, θ ) corresponding to the axonal axis, radial and circumferential directions, respectively, is adopted. The AP and MP coexist in this coupled mechanoelectrophysiological system. The AP is modeled as an electrophysiological pulse using the classic H-H theory [1] (for unmyelinated axons and the Nodes of Ranvier of myelinated axons) combined with the cable theory (for the internodal regions of myelinated axons), whereas the MP is modeled using the viscoelasticity theory in a dynamics framework. Changes of voltage field during an AP induce strain gradient fields on the axon, following the reverse flexoelectric effect [35,49,50], which cause the axon to deform. Reversibly, mechanical deformations in the axon due to a MP not only alter the electrical properties of the membrane (e.g., membrane capacitance and resistance) but also result in a self-polarized electrical field following the direct flexoelectric effect, which in turn influences the membrane potential field. In the following, two types of axons are considered: unmyelinated and myelinated. An idealized representation of these two types of axons with the adopted cylindrical coordinate system is shown in Fig. 1.

B. Electrophysiological pulse model
The electrophysiological pulse model adopted here is an axisymmetric membrane model. This implies a constant voltage throughout the membrane thickness as well as around the axon circumference. Therefore, only the variation of voltage field along the longitudinal direction is considered. The classic H-H model combined with cable theory is used to solve the following partial differential equation describing the propagation of AP along the axon: where t is the time, V is the membrane potential, and I P is the self-polarized current density induced by the membrane deformation due to direct flexoelectricity (see Sec. II D 2). d is the axon diameter which varies with the longitudinal coordinate z. Calculation of d is introduced in Sec. II E. The four terms A, B, C, and D are provided in Table I. The physical definition and values of these parameters are provided in Table III.
In Table I, the dimensionless gating variables m, h, and n follow the evolution equations summarized in Table II, controlled by the six voltage-dependent rate constants α m , β m , α h , β h , α n , β n . See Hodgkin and Huxley [1] for further information.

C. The mechanical pulse model
In the mechanical model, the axon is taken as an isotropic and homogeneous cylinder with a thin wall of constant thickness H. By virtue of the axisymmetry, the displacement field is only a function of the longitudinal coordinate z. The mechanical model starts with the balance of linear momentum:  where ρ is the mass density, v is the velocity, σ is the Cauchy stress tensor, and b is the body force per unit mass. It is observed that the lipid bilayer exhibits both elastic and viscous properties, and its mechanical property may influence the propagation (elastic) or attenuation (viscous) of mechanical signals across the cell membrane [51]. Therefore, we consider two material laws: (1) linear elastic and (2) linear viscoelastic. For (1), σ relates to the infinitesimal strain tensor by: where C is the fourth-order stiffness tensor. For (2), a twobranch generalized Maxwell model is adopted, see Fig. 2. In this case, starting from the general integral representation of linear viscoelasticity, σ relates to through where the internal stress variable h 1 (t ) is defined as: where γ 1 = E 1 E 0 is the normalized elastic modulus for which E 0 and E 1 are the elastic moduli associated to the springs of the purely elastic and viscoelastic branches, respectively, and τ 1 = η 1 E 1 is the relaxation time, where η 1 is the viscosity of the dashpot. These mechanical parameters are calibrated by fitting the equivalent storage and loss moduli of the model to the frequency-dependent rheological properties of neurons experimentally measured in Ref. [52]. Equation (4) is solved in a discretized manner where h 1,n+1 at the t n+1 time step is computed recursively following: where t is the time increment. See Refs. [53,54] for more details on the implementation. The final parameters are provided in Table III.

The reverse flexoelectric effect
To model the reverse flexoelectric effect, it is assumed that a local alteration in the potential of the dielectric membrane from its equilibrium causes a local strain gradient field under 032406-3 the form of either downward bending (membrane depolarization) or upward bending (membrane hyperpolarization) by means of radial force. In addition, it is assumed that the magnitude of membrane bending force is proportional to the offset of the membrane potential from the equilibrium level, i.e., a more polarized membrane leads to a proportionally larger membrane bending force, whereas a less polarized membrane results in a proportionally smaller bending force, see Fig. 3: where F r is the reverse flexoelectric force, f r is the reverse flexoelectric coefficient, and V is the local change in the membrane potential.

The direct flexoelectric effect
To model the direct flexoelectric effect, it is assumed that a strain gradient field in the form of membrane bending can induce a polarization current through the membrane, see spatial gradient of the strain field following: where F is the fourth-order direct flexoelectric tensor and ∵ is the triple-dot product (summing up the two tensors' components on the three inner indices). The polarization current density I P can then be computed as the projection of the rate of change of P on the direction normal to the axon membrane: wheren is the membrane normal. In our study, F reduces to a scalar direct flexoelectric coefficient f d due to the isotropy of the material.
FIG. 4. Illustration of the direct flexoelectric effect: A strain gradient field, e.g., membrane bending, induces a polarization current through the membrane.

E. Alteration of capacitance and resistance due to membrane deformation
The electrophysiological parameters such as membrane electrical constant C m and membrane resistivity ρ m used in Eq. (1) and Table I are fundamental electrical properties of the membrane which only depend on the constitution of the axonal membrane and not its geometry. However, the derived membrane capacitances c m , membrane resistances r m , as well as axial resistances r a used in the finite-element computation do vary with the deformation of the membrane. For an axon with a reference radius of R, the axon diameter d at the location z is calculated by: where w(z) is the radial displacement at the location z and is extracted from the results of the mechanical model. For one element length l of axonal membrane in the deformed configuration, the membrane capacitance c m can be therefore approximated by (assuming constant membrane thickness H m ): Similarly, the membrane resistance r m can be approximated by: The axial resistance r a can be approximated by: In the finite-element implementation, values of the axon diameter d (z) are evaluated at the Gauss points at each time step. Values of c m , r m , and r a are updated for each element accordingly. The same procedure is used for the myelinated membrane.

F. Model coupling
The full coupling of the model is illustrated in Fig. 5. The electrophysiological pulse model is coupled to the MP model via the reverse flexoelectric effect. Conversely, the MP model is coupled to the electrophysiological pulse model through two channels: (1) the direct flexoelectric effect and (2) changes in the electrical properties of the membrane due to membrane deformation, e.g., membrane capacitance and resistance in the H-H model. The coupling parameters are provided in Table III.

G. Parameter selection
The modeling parameters used in our study are summarized in Table III. The values of the geometrical and electrophysiological parameters are taken from the work of Jérusalem et al. [42], in which guinea pig white matter was chosen as the reference material. As the literature data were not directly available for this specific animal and/or brain region, an educated guess based on the similarity with other materials or general scaling laws is made. Original references of the values of these parameters can be found in Jérusalem et al. FIG. 5. Coupling between eletrophysiological and mechanical models. [42]. The values of the mechanical parameters are calibrated by fitting the equivalent storage and loss moduli to the Prony series model of frequency-dependent rheological properties of neurons measured experimentally by Ayala et al. [52]. In addition, the value of the reverse flexoelectric coefficient is determined so that the simulated membrane displacement falls within the range of experimentally observed values, whereas the value of the direct flexoelectric coefficient is determined so that the computed polarization current is of the same order of magnitude as the current source that is used experimentally to initiate an AP. Due to the lack of quantitative information of these parameters in the literature for neuronal membrane, we assume that these values are the same for both Node of Ranvier and internodal regions.

H. Numerical implementation
The axonal membrane is discretized with two-node longitudinal axisymmetric finite elements, representative of the membrane and the underlying cytoplasm in the case of the electrophysiology. Each node is enriched with four degrees of freedom (DoF) corresponding to the voltage, longitudinal and radial displacements, and rotation angle.
The discretized weak form of the electrophysiological governing Eq. (1) can be written as where C el and K el are the electrophysiological damping and stiffness matrices, F el is the electrophysiological external force vector, and V is the unknown nodal voltage DoF vector. The global matrices and vectors are assembled from the element counterparts. A 1D natural coordinate system is 032406-5 adopted within each element so that the domain of each element is transformed to [−1, 1]. The element matrices are evaluated at each time step using four Gauss points and are assembled into the global matrices, see Appendix for the full formulation. The boundary conditions are imposed such that one end of the axon is maintained at the clamped voltage and the other end is modeled as a sealed end where the spatial longitudinal gradient of V is 0. Solutions at the nodal points are obtained using the explicit forward difference scheme. Finally, values of m, h, and n are updated as internal variables stored at the Gauss points using the corresponding equations in Table II.
The discretized weak form of the mechanical governing Eq. (2) can be written as for a linear elastic material, where M m and K m are the mechanical mass and stiffness matrices, F m is the mechanical external force vector, and d is the unknown nodal displacement DoF vector. Note that d is defined for a 1D two-node shell element, and each node is assigned a longitudinal displacement, a radial displacement, and rotation angle. A 2D natural coordinate system (ξ, η) is adopted. In this work a linear viscoelastic material model is adopted and Eq. (16) needs to be evaluated iteratively between time step t n and t n+1 : where K vis m and K hist m are the viscous stiffness and history stiffness matrices, and H hist m,n+1 is the history matrix at the t n+1 time step. Note that the evaluation of element matrices in the two equations above requires volume integrations over a differential volume. However, because of the structural axisymmetry, the integrands of the element matrices are all independent of the circumferential coordinate θ , and the volume integrations can be converted to double integrations. These element matrices are evaluated at 2×4 Gauss points at each time step and are assembled into the global matrices. For the viscoelastic material model, evaluation of the element history matrix at the t n+1 time step requires the value of the internal variable h 1,n at the previous time step, and its value is recursively updated following Eq. (7). After assembly of the element matrices, boundary conditions are applied as a pinned joint at one end and as a roller support at the other end. The force vectors include the contribution of the direct flexoelectric force following Eq. (8), see the Appendix for more details. The Newmark method with values β = 0.25 and γ = 0.5 is used as the time marching scheme to obtain the solution at each time step. The updated displacement (strain) field at each time step is used in (a) Eq. (1) to update the axon radius R, (b) Eq. (9) to compute the polarization density, and (c) Eqs. (7) and (17) to update the terms which are related to the stress and strain history in the viscoelastic material model.

III. RESULTS
Numerical simulations are conducted for several different scenarios: (a) comparison of the difference of AP conduction between unmyelinated and myelinated axons, b) comparison of MP below and above "activation threshold", (c) pulse collision between two electrical pulses (EPs) 1 as well as between an EP and a MP, (d) pulse chasing between an EP and a MP, and (e) effect of changing the pulse repetition frequency (PRF) of the mechanical wave. Each of the results is illustrated and discussed separately in the following subsections. In all following simulations where a MP is triggered, its period is 0.2 ms.
A. Unmyelinated vs. myelinated Figure 6 illustrates the propagation of two EPs both initiated by electrical stimuli on an unmyelinated axon and on a myelinated axon. The corresponding video is provided as Supplemental Material [55]. In both cases, an accompanying MP can be observed to travel conjointly with the EP, with a magnitude of approximately 1 nm. This MP is linked to the axon surface movement occurring during membrane polarization due to the reverse flexoelectric effect. In the unmyelinated axon, the axon is fully covered by ion channel regions, which, when activated by a local elevation of voltage, activate their neighboring regions in the same way. As a result, all ion channels are activated in a continuous sequence and the EP travels forward maintaining a wave shape. On the contrary, in a myelinated axon, the axon is covered intermittently by regions of Nodes of Ranvier and internodal regions. When one nodal region is activated by a locally elevated voltage, the membrane potential in the neighboring internodal region slowly propagates to the next node in a quasidiffusion fashion. When the voltage at the end of the internodal region is elevated enough, the next nodal region is activated. As a result, the electrical potential thus propagates segment by segment as shown in Fig. 6 rather than in the continuous wave form as in an unmyelinated axon. Figure 7 shows the propagation of two MPs in a myelinated axon, both of which are initiated by stimulating the right end of the axon with a waveform displacement profile. The corresponding video is provided as Supplemental Material [55]. The two cases differ through the magnitude of the stimulating displacement profile: One involves a sudden overall compressive strain of 0.5% and the other of 0.25%. In the first case, the polarization current induced by the local strain gradient due to propagating mechanical wave is large enough to activate the ion channels' activities: An EP is initiated and propagates conjointly with the MP, i.e., at the velocity of the MP as opposed to the velocity of a sole EP. In the second case, however, the induced polarization current is too small to activate the ion channels' activities. As a result, only the propagation of a local small membrane polarization could be observed. This result indicates the existence of an "activation threshold" which corresponds to the smallest magnitude of the mechanical stimulus that is needed to activate the ion channels' activities through the flexoelectric effect. FIG. 6. Comparison between the conduction of an AP on an unmyelinated axon vs. on a myelinated axon. In both cases, an accompanying MP can be observed to propagate conjointly with the EP with a magnitude of approximately 1 nm. The corresponding video is provided as Video 1 in Supplemental Material [55].

B. Above vs. below threshold
FIG. 7. Comparison between the conduction of a MP with magnitude above or below the "activation threshold" on a myelinated axon. In the case of MP above the threshold, the MP is able to induce an EP due to the direct flexoelectric effect which conjointly propagates with the MP. On the contrary, in the case of MP below threshold, the flexoelectric current induced by the MP is too small to trigger an EP. The corresponding video is provided as Video 2 in Supplemental Material [55].
FIG. 8. Collision between an EP and a MP on a myelinated axon. The left end is excited by an electrical stimulus, whereas the right end is excited by a mechanical stimulus with 9-ms time delay. The magnitude of the MP is large enough to induce a second EP. It is observed that the two APs annihilate each other up on collision, whereas the MP maintains itself and induces a new EP after leaving the refractory period zone of the right-moving EP. The corresponding video is provided as Video 3 in Supplemental Material [55]. Figure 8 shows the collision between an EP and a MP in a myelinated axon. The corresponding video is provided as Supplemental Material [55]. The EP is initiated by an electrical stimulus on the left end of the axon, whereas the MP is initiated by a mechanical stimulus on the right end whose magnitude is above the "activation threshold." The leftward propagating MP itself does not interfere with the rightward propagating EP. Nevertheless, the mechanically induced leftward propagating EP interferes with the original rightward propagating electrical pulse. After collision, these two EPs annihilate each other, and all electrical activities temporarily die out, whereas the MP continues to propagate. After leaving the refractory period zone of the rightward propagating EP, the MP is free to trigger another EP, which can copropagate with the MP. If both the simulation time and length of the axon are long enough, then the MP is gradually dissipated by the viscosity of the axon membrane and eventually vanishes. Another simulation (results not shown) was done on the collision of two MPs, both of which able to generate a flexoelectric current strong enough to induce an EP. In this situation, it is similarly observed that the two mechanically induced EPs annihilate after collision, whereas the two MPs themselves pass each other without any mutual influence. These two MPs eventually vanish due to the viscous dissipation of their kinetic energy. Figure 9 shows an EP chased by a MP in a myelinated axon. The corresponding video is provided as Supplemental Material [55]. To this end, an EP is initiated from the right FIG. 9. Chasing of an EP by a MP. The EP is initiated by an electrical stimulus from the right end. After a 2-ms time delay, the MP is triggered from the right end with a magnitude above the "activation threshold" to induce a second EP. The MP carries the second EP to propagate at the speed of the MP. However, when the second MP enters the refractory period of the first EP, it is split from the MP and gradually dies out, whereas the MP continues to propagate, catches up with the first EP, pushes it to propagate at the speed of the MP, and enhances the magnitude of the EP. The corresponding video is provided as Video 4 in Supplemental Material [55]. For case 1, the conduction of the EP on an intact "healthy" axon triggered by an electrical stimulus on the right end is used for comparison. For cases 2-6, a Node of Ranvier is intentionally myelinated in an idealized case of damaged node. In case 2, the axon is excited by an electrical stimulus from the right end; the EP cannot pass the damaged node and eventually vanishes. In case 3, a single MP (smaller than the "activation threshold") is sent 2 ms after the EP. While the MP itself can propagate through the damaged node, it is unable to carry the EP across it. In case 4, MPs are triggered at 1 kHz, and this series of MPs helps the EP to conduct further away, but the EP still gradually diminishes and is not able to activate the next Node of Ranvier. In case 5, the PRF of the MPs is increased to 5 kHz, and now the MP is able to carry the EP over the damaged node. However, the conduction velocity of the EP in this case appears as reduced by the damaged node in comparison with case 1 of the "healthy" axon. Finally, in case 6, the PRF of the MPs is increased to 10 kHz. It is observed that the MPs not only carry the EP over the damaged node but also enhance both the polarization level as well as the conduction velocity of the EP. The corresponding video is provided as Video 5 in Supplemental Material [55]. All snapshots are at time 7 ms. end by an electrical stimulus. After 2 ms, a MP is also initiated from the right end whose magnitude is above the "activation threshold" and therefore capable of inducing a second EP. This second EP is carried by the MP propagating at a speed faster than the first EP. However, when arriving at the next Node of Ranvier, the second EP is split from the MP as the node is still within the refractory period of the first electrical pulse. From here, the MP propagates at its own speed, whereas the second EP gradually dies out. When the MP catches up to the first EP, it accelerates the first EP to propagate at the speed of the MP but also enhances the level of membrane polarization.

E. Frequency-dependent neuromodulation below activation threshold
In this study, the effect of changing the PRF of the MP is also briefly investigated in the context of electrophysiological enhancement of damaged myelinated axons in Fig. 10. For comparison, simulation of conduction of an EP on an intact "healthy" axon is also included (case 1). To simulate a damaged myelinated axon, the dynamics of one Node of Ranvier is stopped by integrating the corresponding elements into their two neighboring internodal regions. The cable theory used to model the internodal regions is naturally associated with a length constant (i.e., the length of the internodal region), which is representative of the distance that a graded electric potential can travel along the axon through passive electrical conduction while still being able to activate the next Node of Ranvier. As a result, when the internodal region is too long (as in this case, where two consecutive internodal regions are concatenated), the membrane potential is unable to travel through, and the next Node of Ranvier is not activated (case 2). For this scenario, the simulation of Sec. III D is first repeated, but with a MP whose magnitude is lower than the "activation threshold" (case 3). In this case, the MP is unable to carry the electrical potential to "jump" over the damaged Node of Ranvier region and the electrical potential gradually dies out. Another similar simulation is also conducted with a series of MPs at 1 kHz PRF and an "infra-activation threshold" magnitude (case 4). It is shown that this series of MPs helps the EP to conduct further away, but the EP still gradually diminishes and is not able to activate the next Node of Ranvier. However, if the PRF of the MPs is increased to a level above 5 kHz (cases 5 and 6), then it is observed that the MP train under this frequency not only can help the EP to jump over the damaged region but also may enhance its propagation in the context of both polarization level and conduction velocity, even if the magnitudes of each MP taken individually are all below their individual "activation threshold," see Fig. 10 and Supplemental Material [55].

IV. DISCUSSION
Mechanical deformation of the membrane during membrane electrical activity has been experimentally observed in a wide variety of excitable cells and tissues including both invertebrate and vertebrate nerve fibers [2][3][4][6][7][8][9][10][11][12][13][14][15][16][17][18][19]21,22]. The magnitude of the displacement varies depending on the type of cell and nerve fiber, e.g., about 0.5-1 nm in squid giant axons, 5-10 nm in crab nerves, and 10-20 nm in crustacean nerves. In addition, the magnitude also appears to be proportional to the change of voltage, in the ratio of about 1 nm per 100 mV. In recent years, numerical models have been proposed to study this phenomenon. Among others, Hady and Machta [43] proposed a numerical model in which these copropagating mechanical displacements (termed "action waves") emerge from the surface waves due to the varying compressive electrostatic forces across the membrane, whereas these varying electrostatic forces are induced by the traveling wave of the electrical depolarization. Their simulation results show a good agreement with a range of experimental results on the magnitude of the membrane displacement and the wave speed. Interestingly, they concluded that many properties of the comoving mechanical wave are determined by comparison of its own propagation velocity and the velocity of the EP at relevant wavelengths. Drapaca [44] also proposed a coupled electromechanical model of a neuron based on a unified variational principle. In this model, the neuron is modeled as a linear viscoelastic Kelvin-Voigt solid whose electrochemical activity is described by the classical H-H equations. The neuron electrodynamics is described by a coupled system of differential equations, which is obtained by minimizing a special integral functional whose integrand is made of the kinetic and potential energies as well as the work done by the forces acting on the neuron. The membrane capacitance is also linked to the mechanical deformation of the neuron. However, the elastic modulus of the neuron is not constant. Instead, it is assumed to depend on the gating variables in the H-H model. Finally, the AP is simulated by applying a constant external electric current and an initial displacement as well as an initial velocity of the membrane. In a series of publications [26][27][28], Engelbrecht et al. extensively studied the electromechanical coupling of waves in nerve fibers. In their model, the AP is modeled as an EP which is assumed to trigger all other processes. It is modeled using a modified FitzHugh-Nagumo model [56,57] which also incorporates "mechanical activation" parameters that describe the mechanosensitivity of the ion channels. The EP is coupled to the mechanical waves which include both longitudinal and transverse waves in the membrane, as well as the pressure wave in the axoplasm. The longitudinal wave is modeled using the improved nonlinear soliton model, whereas the transverse wave is linked to the longitudinal wave via the theory of rods. The pressure wave in the axoplasm is modeled using a model of pressure propagation in an elastic cylindrical tube with added viscous damping term described by a 1D Navier-Stokes equation. The coupling effect is achieved via several coupling forces added into the governing equations. These coupling forces depend on the spatial and temporal changes of voltage. Apart from these models, there are also models which attempt to explain the membrane polarization-induced deformation based on the equilibrium of cell membrane tension, e.g., cells alter their radii to attempt to maintain constant pressure across the membrane due to the change in surface tension [10,16,58]. In these references, the surface tension is defined by the Gibbs free energy required to maintain a certain surface area, containing both chemical and electrical contributions. The chemical contribution is the natural surface tension of membrane interface with no electric field, whereas the electrical contribution comes from the electrical energy stored in the membrane through the membrane capacitance during polarization induced by the applied voltage field. When the membrane is polarized, the change of membrane potential modulates the surface tension of both intracellular and extracellular fluids with the deformable cell membrane. In order to maintain a constant pressure across the membrane, differences in surface tensions between the intracellular and extracellular interfaces induce changes in membrane curvature, exhibited as membrane displacement. Consequently, changes in membrane potential result in an alteration of cell radii in the form of membrane movement to maintain the mechanical equilibrium across the membrane [59]. In the study presented here, the EP is modeled using the class H-H theory for the ion channel mechanism with the cable equation for its conduction, whereas the MP is modeled using the dynamic beam assumption for a viscoelastic solid. The effect of the pressure wave in the axoplasm is ignored. The membrane polarization-induced deformation is captured by the reverse flexoelectric effect, where a local potential change always induces a membrane displacement in the form of a local strain gradient field. This membrane deformation can propagate independently as a mechanical wave, whose speed depends on the mechanical properties of the membrane. If the membrane is stiff enough, then it may actually propagate ahead of the EP. However, because of its low magnitude and the viscoelasticity of the membrane, it is eventually damped and vanishes. Simulation results of our study also show good agreement with the experimental measurements, i.e., about 1-nm membrane displacement.
The conduction velocity (CV) of an AP in the axons can be as slow as 0.07 m/s for specific frog muscle fibers [60] and as fast as 120 m/s for Type Ia sensory fibers [61], with typical values ranging from 1 m/s to over 100 m/s [62]. The CVs are known to be affected by many factors, such as species, gender, age, type of nerve, axonal diameter, and, most importantly, axonal myelination. Myelin sheaths on the neuronal axons can help the nerve system to acquire a faster propagation of the AP while increasing energy efficiency by means of saltatory conduction [63]. Normally, the conduction velocity of an AP in the myelinated axons can be several folds faster than that in the unmyelinated axons [64]. In the simulations presented here, the CVs of EPs without MPs (as MPs drive EPs' velocities) are approximately 0.63 m/s in both myelinated and unmyelinated axons. This value falls in the typical value range of CV of an AP. However, the model does not predict a larger CV in myelinated axons when compared to the unmyelinated case. This could be attributed to the following reasons. First, in the proposed modeling framework, the classic cable theory (initially developed by Thomson [65] to model signal decay in submarine telegraphic cables) is used to model the propagation of electrical signal along the passive myelinated axon. It considers the axon as a cylinder composed of individual segments of membrane capacitances and resistances aligned in parallel, intermittently connected by the axoplasm resistances, and the myelin layers as a series of myelin capacitances in serial with the membrane capacitor, as well as a series of resistances in serial with the membrane resistances. This type of model has been used widely in the literature to model the propagation of an EP in a myelinated axon. However, given the complex dual nature of AP discussed in the paper, the actual physiological mechanisms by which myelin layers accelerate the EP propagation cannot be reduced to a combination of simple electronic components. Second, as stated in Sec. II G, the guinea pig white matter was chosen as the reference material for evaluating the electrophysiological parameters. As the literature data of the membrane and myelin electrophysiological parameters for this specific type of axons are not directly available, educated guesses based on the similarity with other materials from different references along with general scaling laws are made. Therefore, using values of electrophysiological parameters from a combination of several sources may not necessarily produce a result that can perfectly simulate the conduction of an EP on a myelinated axon in every aspects. Finally, since the CV of an AP in myelinated axons also vary between different cells and nerves, potentially depending on the type and population of ion channels, axon diameter, thickness, and maturation of the myelination, ideally each cell or axon should be submitted to a thorough histological study identifying the population of ion channels and myelination properties. Based on this, the electrophysiological parameters tailored to this specific cell or axon type could be calibrated. However, this type of work itself is a full enterprise on its own (e.g., one of the main aims of the Blue Brain Project is to gather ion channel models, see http://channelpedia.epfl.ch/), but, while acknowledging the limitations of these models, the qualitative conclusions drawn from these simulations should not be a priori affected by the CV of the EPs. Figure 8 shows the collision between an EP and a MP in a myelinated axon. During the simulation, the left-moving MP induces a left-moving accompanying EP which leaves a refractory period behind. When the right-moving EP collides with this mechanically induced left-moving EP, it enters the refractory period of the left-moving EP. Since the ion channels in the refractory period zone cannot be immediately activated again, the right-moving EP cannot propagate through and therefore dies out. However, the MP itself is free to propagate. After leaving the refractory period zone of the right-moving EP, the MP again induces polarization in the membrane due to the direct flexoelectric effect. If the viscosity of the membrane is not set to be large enough to have attenuated the MP below its "activation threshold" (as shown in this case), the membrane polarization can develop into a new EP which copropagates with the MP. This phenomenon appears even more clearly in the case of unmyelinated axons, in which the postannihilation re-induced membrane polarization may develop and split into two EPs separating and propagating in opposite directions.
In the literature, several experimental studies were conducted in various nerve fibers or cells to examine whether annihilation can occur or not on AP collision [66][67][68][69][70]. It has been frequently witnessed that colliding APs are reciprocally annihilated instead of penetrating each other. One exception is a series of studies by Gonzalez-Perez et al. [18,19], where it is observed that the two simultaneously generated pulses propagating in orthodromic and antidromic direction in the giant axons of both earthworms and lobsters pass through rather than annihilate each other. The numerical simulation results illustrated in these studies where the AP is modeled as a mechanical wave based on the soliton theory also support their experimental observations. Nevertheless, in one comment to the original paper [71], the experiment was repeated but the results about AP penetration could not be reproduced as AP annihilation on collision was consistently observed. In the following reply to this comment [72], it was argued that the signals had been strongly perturbed due to the significantly different experiment configuration in the study of Berg et al. [71]. In our study, it is observed that when two EPs from opposite directions meet and collide, they annihilate each other regardless of whether the EPs are initiated by electrical stimuli or triggered by MPs whose amplitudes are above the "activation threshold." These results are consistent with most of the experimental observations on nerve pulse collision. The fundamental reason is specifically attributed to the mathematical model of the EP that is adopted in our model, i.e., the H-H model. Due to the concept of refractory period in the H-H model, the ion channels are relaxed and, therefore, cannot be activated again within a short time. Therefore, when two EPs collide, they both enter each other's refractory period, hence their mutual annihilation. The soliton hypothesis, on the contrary, considers APs as soundlike MPs and predicts no annihilation of colliding pulses (as shown in the studies of Gonzalez-Perez et al. [18,19]). However, most recently, Shrivastava et al. [73] observed the annihilation of two superthreshold interfacial MPs on head-on collision in a lipid monolayer. The blatant discrepancy between the prediction of theoretical analysis of the soliton hypothesis and the experimental observations of MP collision was attributed to the result of nonlinear material properties (such as nonlinear viscosity resulting from relaxation of phase change), as well as the internal heat transfer in the colliding pulses. Consequently, it is clear that the most critical factor is still the predefined conception toward the nature of an AP. In our study, we assume that APs are merely electrophysiological waves due to the ion channel mechanism, whereas the associated mechanical wave corresponding to the membrane displacement is a consequence of this electrical phenomenon due to flexoelectricity. This conception certainly leads to simulation results in which wave annihilation is expected on collision. To unveil the real nature of APs, experimental measurements of MPs and their interactions in the lipid bilayers of cell or nerve membranes (if technology permits) is urgently needed. However, the approach we are taking essentially demonstrates that the membrane displacement, observed by the proponents of the soliton theory, is not necessarily incompatible with the H-H model when one accounts for the flexoelectricity of the membrane.
Many sophisticated models of biological flexoelectricity have been proposed in the literature for the purpose of electromechanical coupling. For example, Rey [74] developed a formulation of an isotropic interfacial liquid crystal flexoelectric membrane under different mechanical deformation conditions and derived a membrane electromechanical shape equation that connects fluid forces with membrane curvature and electric displacement. This model is incorporated and further sophisticated by Rey et al. [75] to model a mechanical energy harvesting system consisting of a deformable soft flexoelectric thin membrane, as well as by Herrera-Valencia and Rey [76] to model the actuation of flexoelectric membranes of hair cells in viscoelastic fluids. Gao et al. [38] developed an electromechanical liquid crystal model for characterizing the mechanical equilibrium morphology of an axisymmetric lipid vesicle in a uniform electric field. In this model a general equation was established to govern the vesicle shape based on the Helmholtz free energy, which incorporates the effects of elastic bending, osmotic pressure, surface tension, Maxwell pressure, as well as flexoelectric and dielectric properties of the lipid membrane. Breneman et al. [77] proposed a flexoelectric model for stereocilia of the hair cells, in which the governing equations include the electrical cable and mechanical wave equations in which the corresponding terms are replaced by their flexoelectric counterparts (i.e., flexoelectric current and flexoelectric axial stress). Mohammadi et al. [78] proposed a theory for flexoelectric membranes based on minimization of the total free energy as the sum of the internal and electric energies, both of which are functions of out-of-plane displacement and out-ofplane polarization, for a 2D thin membrane. This model was extended by Deng et al. [40] to model flexoelectricity in soft materials. More recently, Kancharala et al. [79] developed a flexoelectric model coupled with a nonlinear finite-element model to study the dynamic mechanoelectrical response of droplets under excitation. This electrostatic model captured the internal distribution of charge within a droplet, leading to a dipole and surface potential. In our study, two simple phenomenological models are used to model the reverse and direct flexoelectric effects, see Eqs (8) and (9). Values of the reverse and direct flexoelectric coefficients are calibrated so that the simulated membrane displacement falls in the range of experimentally observed values and the computed polarization current is of the same order of magnitude as the current source that is used experimentally to initiate an AP. The electrophysiological and mechanical models are mutually coupled via these two simple phenomenological flexoelectric models. While the approach would benefit from more complex considerations accounting for the fundamental mechanisms that underly flexoelectricity [29,33,35,40,49,50,80,81], it provides a flexible tool for the investigation of the coupling between mechanics and electrophysiology in the context of neuromodulation.
In recent years, the potential of low-intensity, lowfrequency transcranial pulsed ultrasound (TPU) to stimulate the brain as a noninvasive neuromodulation method has drawn the interest of many research groups. Compared with more traditional neuromodulation methods such as deep brain stimulation or vagal nerve stimulation, TPU does not require the implantation of electrodes that could damage the nervous tissue or lead to infection. In many studies, TPU has been shown to induce changes in the animals' activities and neural response [82][83][84][85][86][87][88][89][90][91]. Moreover, TPU has also been experimentally observed to affect the neural activity of human beings [92][93][94] by suppression or inhibition dependent on the parameters of the acoustic energy applied to neural tissue. Systematic reviews of TPU for neuromodulation can be found in Refs. [95][96][97][98][99][100][101][102][103]. Despite these experimental evidences, remarkably, the action mechanism underlying ultrasound neuromodulation is still poorly understood.
In an attempt to answer this question, several hypotheses have been proposed. For instance, Tyler [104] suggested that ultrasound alters the membrane conductance by interacting directly with the neuron membrane and its surrounding fluid environment via various types of mechanical actions. Plaksin et al. [105] proposed a model in which ultrasound-induced intramembrane cavitation within the bilayer membrane leads to cell excitation through the effect of currents induced by membrane capacitance changes. In this study, an idealized loading case is considered: A temporally varying strain is applied at one axonal end as a boundary condition, resulting in a traveling local strain gradient. For any given frequency, the applied strain level is incrementally increased to identify an "activation threshold" above which a single or multiple MPs initiate an EP by direct flexoelectric effect. The simulation results demonstrate that such MP-induced EP can potentially either suppress or enhance an EP, depending on the directionality of the two pulses. In classical TPU neuromodulation, the acoustic wave is usually frequency, intensity, and power modulated, and the propagating wave in the skull is attenuated due to the scattering and absorption of the brain tissue [106]. Therefore, there is no direct experimental observation of the level of physiological deformation experienced by a neuron during TPU neuromodulation. As such, while both idealized and decoupled from what would potentially be required at the tissue level, the proposed loading scenario has the merit to focus on the neuromodulation at the cell scale. It must also be emphasized that, according to the direct flexoelectric effect, it is the local strain gradient field over a very short moment that leads to the initiation of an EP, not the strain or displacement fields [see Eqs. (9) and (10)]. Therefore, in this theory, normal body movement or local musculature changes are not expected to trigger an AP, as either the neurons do not experience a gradient of deformation (such as for homogeneous deformation or displacement due to head movement) or the duration under which the neurons experience the strain gradient field is too long to trigger an EP.
The simulation results also show that, when emitted at a given PRF, MPs with small magnitudes (even below the single MP "activation threshold") can also enhance the AP and help it propagate through a damaged region of axon. The natural frequency of these MPs is set to be 5 kHz which is much lower than the natural frequency of a TPU. One additional simulation (see Video 6 in Supplemental Material [55]) of mechanical-induced EP was conducted on a "healthy" myelinated axon stimulated by a train of MPs with a natural frequency in the range typically used in experimental TPU [107]: 500 kHz for 1 ms for a full simulation of 3 ms, i.e., a duty factor of 33.3%. The amplitude of the MPs used in this case is much lower than that used in the low frequency cases (approximately 1/10th). However, the train of MPs at this natural frequency is capable of inducing an EP which includes both depolarization and repolarization phases. This result indicates that while a single wave requires a relatively large displacement (approximately 1 μm), much lower membrane displacements are required for neuromodulation at high frequency. It can also be seen that this train of MPs can induce an EP larger in amplitude when compared with the ones in Fig. 10. Therefore, a similar result can be expected: A MP with a typical TPU frequency can also enhance the EP and help it to propagate through the damaged region on the axon. To fully examine the correlation among PRF, natural frequency, and amplitude, additional simulations are needed, which simulate EPs triggered by MPs for different frequencies and amplitudes. While it is out of the scope of this study, it will be investigated in the near future.

V. CONCLUSION
In order to investigate some of the modern multiphysics hypotheses on the nature of AP, we proposed a mechanoelectrophysiological coupling thin-wall model for axons in neurons. The electrophysiological component of the AP is modeled as an electrophysiological pulse whose conduction is described by the cable theory in the myelinated regions, and the H-H equations at the Nodes of Ranvier and for unmyelinated axons. The mechanical response of the axon is modeled using the beam theory for a viscoelastic material applied to an axisymmetric geometry in a dynamic framework. The electrophysiological and mechanical models are coupled via the direct and reverse flexoelectric effects. The deformation of the membrane also alters the membrane electrical properties through the change of the modeling parameters in the H-H model. Using this model, a series of simulations was systematically conducted on both myelinated and unmyelinated axons. It is shown that an AP is always accompanied by an in-phase propagating membrane deformation; and a MP with enough magnitude can also initiate an AP due to the flexoelectric effect. The results of pulses interference depend on their respective directionalities. Independently increasing the amplitude or PRF of the MP can enhance the propagation of an AP, even on a damaged axon. Our study provides a potential explanation to these experimentally observed axon membrane displacements during AP propagation. It also provides a foundational framework for the modeling of the mechanisms behind ultrasonic neuromodulation.
The code used for the finite-element formulation is available in Ref. [108].

Finite-element formulation of the electrophysiological model
The strong form of the electrophysiological governing equation reads: Its finite-element discretized weak form at the element level can be written as: where two-node line elements are adopted with V e = [V i V j ] T being the element degree of freedom at the two end nodes i and j. A 1D natural coordinate ξ with −1 ξ 1 is adopted for each element. The local coordinate z within each element is mapped to the natural coordinate ξ using the coordinate transformation: z = l 2 + l 2 ξ with the Jacobian J el = ∂z ∂ξ = l 2 and its inverse J −1 el = ∂ξ ∂z = 2 l , where l is the length of the element. The element shape function matrix N el is given by: The derivative of N el with respect to the natural coordinate ξ, B el , is then given by: The element damping matrix C e el is given by: while the element stiffness matrix K e el is given by: and the element force vector F e el is given by: where Within each element, values of C e el , K e el are evaluated using Gaussian quadrature using four Gauss points at each time step: where ξ i denotes the ith Gauss point's natural coordinate and W ξ i is its corresponding weight. Finally, the element force vector F e el are evaluated by: There is no need to compute the actual values of Q ep , as when the element force vectors are assembled into the global force vector F el , the values of Q ep in the two adjacent elements cancel each other. The explicit forward difference scheme is used here. The initial condition is that the membrane potential at time t = 0 is equal to the membrane resting potential V (t = 0) = V rest . The boundary conditions differ between the different modeling scenarios. At each time step, the element matrices C e el and K e el as well as the force vector F e el are assembled into global matrices C el , K el and global force vector F el , respectively: The solutions at the nodal points are obtained and interpolated at the Gauss points using the shape functions to update values of m, n, and h for the next time step.

Finite-element formulation of the mechanical model
The strong form of the mechanical governing equation reads: where two-node axisymmetric shell elements are used with d e = [u i w i θ i u j w j θ j ] T being the element degree of freedom at the two end nodes i and j, where u denotes the axial displacement, w denotes the radial displacement, and θ denotes the rotation angle. We adopt a 2D natural coordinate system (ξ, η) within each element, as shown in Figs. 11 and 12, with ξ being parallel to z and η being parallel to r, and −1 (ξ, η) 1. The local coordinates (z, θ ) within each element are mapped onto the natural coordinates (ξ, η) using the coordinate transformation: z = l 2 + l 2 ξ and r = R + H 2 η, where R is the initial radius of the axon and H is the assumed constant thickness of the membrane. Therefore, the Jacobian matrix J m of the transformation reads: with the Jacobian J m = |J m | = Hl 4 . The element displacement shape function matrix N m is given by: FIG. 11. Axisymmetric differential volume of the axonal cylindrical structure.
with the shape functions expressed in the local natural coordinates: Due to the axisymmetric deformation, the strain tensor only has three nonzero entries (γ θz = γ rθ = 0 and r = 0 as (A16) In axisymmetric cylindrical coordinates, the straindisplacement relationships read: The strain can then be expressed in terms of the nodal displacement field by = B m · d e , where the matrix B m is defined as: The stress-strain relationship for an axisymmetric isotropic linear two-branch generalized Maxwell material reads: where γ 1 = E 1 E 0 is the normalized elastic modulus, where E 0 and E 1 are the elastic moduli of the corresponding spring components. τ 1 = η 1 E 1 is the relaxation time, where η 1 is the viscosity of the dashpot component. D is given by: (A20) Equation (A19) can be formulated in a time-discretized manner as: where t = t n+1 − t n is the time increment and A 1 is defined as: (A22) The internal stress variable h 1 at time t n is computed recursively by: As shown in Fig. 11, a differential volume dv in a cylindrical structure can be expressed as dv = rdθ drdz = J m (R + H 2 η)dθ dξ dη in the natural coordinate. This relationship is used to convert the volume integrals in the following computation of the stiffness and mass matrices to double integrals using the local natural coordinates (as the integrands are independent of θ , the integral with respect to θ is 2π ).
Making use of the equations above, the element mass matrix M e m reads: The element viscous stiffness matrix K vis,e m is given by: The element history stiffness matrix K hist,e m is given by: where ξ i and η j denote the ith Gauss point coordinate in the ξ direction and the jth Gauss point's coordinate in the η direction, respectively. W ξ i and W η j are the corresponding Gauss weights of the ξ i and η j Gauss points. When assembled, the global force vector F m is straightforwardly obtained (without need for assembly) from: where F r n is the reverse flexoelectric force at the nth nodal point. The calculation of F r n is described by Eq. (8).
The Newmark method is used as the time marching scheme, with values of β = 0.25 and γ = 0.5. The initial conditions are set so that the axon is in a static state, i.e., d 0 = 0 andḋ 0 = 0. Note, however, thatd 0 = 0 as the membrane is subjected to initial reverse flexoelectric forces resulting from the membrane potential changes at initial time t 0 = 0. The boundary conditions vary between the different modeling scenarios, depending on which end is stimulated mechanically.
The solutions are obtained at the nodal points and interpolated at the Gauss points using the shape functions to update the internal variables and calculate the stress.