Energy-Optimal Electrical-Stimulation Pulses Shaped by the Least-Action Principle

Electrical stimulation (ES) devices interact with excitable neural tissue toward eliciting action potentials (AP’s) by specific current patterns. Low-energy ES prevents tissue damage and loss of specificity. Hence to identify optimal stimulation-current waveforms is a relevant problem, whose solution may have significant impact on the related medical (e.g. minimized side-effects) and engineering (e.g. maximized battery-life) efficiency. This has typically been addressed by simulation (of a given excitable-tissue model) and iterative numerical optimization with hard discontinuous constraints - e.g. AP’s are all-or-none phenomena. Such approach is computationally expensive, while the solution is uncertain - e.g. may converge to local-only energy-minima and be model-specific. We exploit the Least-Action Principle (LAP). First, we derive in closed form the general template of the membrane-potential’s temporal trajectory, which minimizes the ES energy integral over time and over any space-clamp ionic current model. From the given model we then obtain the specific energy-efficient current waveform, which is demonstrated to be globally optimal. The solution is model-independent by construction. We illustrate the approach by a broad set of example situations with some of the most popular ionic current models from the literature. The proposed approach may result in the significant improvement of solution efficiency: cumbersome and uncertain iteration is replaced by a single quadrature of a system of ordinary differential equations. The approach is further validated by enabling a general comparison to the conventional simulation and optimization results from the literature, including one of our own, based on finite-horizon optimal control. Applying the LAP also resulted in a number of general ES optimality principles. One such succinct observation is that ES with long pulse durations is much more sensitive to the pulse’s shape whereas a rectangular pulse is most frequently optimal for short pulse durations.


Introduction
Electrical stimulation (ES) today is an industry worth in excess of 3 G$. ES devices interact with living tissues toward repairing, restoring or substituting normal sensory or motor function [1]. The rehabilitation-engineering applications scope is constantly growing: from intelligent limb prosthetics and deep-brain stimulation (DBS) to bi-directional brain-machine interfaces (BMI), which are no longer just about recording brain activity, but have also recently used ES toward closed-loop systems, [2][3][4][5].
Application-specific current patterns need to be injected toward reliably eliciting action potentials (AP's) in target excitable neural tissue. To prevent tissue damage or loss of functional specificity, the employed current waveforms need to be efficient. This may significantly impact the biomedical effects and engineering feasibility. Hence, an optimization problem of high relevance to the design of viable ES devices is to minimize the energy required by the stimulation waveforms, while maintaining their capacity for AP triggering toward achieving the targeted functional effects.
A number of recent studies of ES optimality are based on extensive model simulation and related numerical methods through the wider spread of high-performance computing, e.g. [6][7][8][9]. The model dynamics to iterate can be arbitrarily complex and nonlinear. This implies lengthy numerically-intensive computation, irregular convergence and constraints that may be difficult to enforce -e.g. that an AP is an all-or-none phenomenon. Thus, any function of membrane voltage will suffer dramatic discontinuities at parameter-space manifold boundaries where intermittent AP's are likely to be elicited.
Hence, such an iterative approach is not only computationally expensive, but its solution quality is highly uncertain and modelspecific. The long-lasting iteration may converge to shallow local energy-minima. Such numerical misdemeanor of the approach is well known to its frequent users.
In this work we follow the ES pioneers -we use physical reasoning and related mathematics toward a more theoretical treatment of the subject.
Below we summarize very briefly our historical premises. ES' theoretical cornerstones were laid a century ago by experimentally-driven assumptions and models, [10][11][12]. Various constant ES current levels and durations were tried systematically. E.g. Louis and Marcelle Lapicque spent many years performing such lab experiments with multiple physiological preparations [13,14]. This classical work led to concepts like strength-duration curve (SD), i.e. the function of threshold (but still AP-evoking) ES current strength on duration. The first mathematical fit to this empirical results is usually attributed to Weiss, [10,15] where T is the stimulus duration, b is called the rheobase (or rheobasic current level) and c is the chronaxie.
The most expedite way of introducing the rheobase and chronaxie would be to point to eqn. (1) and notice that: i.e. the rheobase is the threshold current strength with very long duration, and chronaxie is the duration with twice the rheobasic current level. In the pioneering studies electrical stimulation was done with extracellular electrodes. Eqn. (1) is the most simplistic of the 2 'simple' mathematical descriptors of the dependence of current strength on duration, and leads to Weiss' linear charge-transfer progression with T, Q(T)~T|I THR~b |(Tzc): Both Lapicque's own writings - [11][12][13], and more recent work are at odds with the linear-charge approximation. Already in 1907 Lapicque was using a linear firstorder approximation of the cell membrane, modeled as a single-RC equivalent circuit with fixed threshold: with time constant t~C=g; C and g~1=R are the membrane capacity and conductance respectively. The second form of eqn. (4) is easily obtained by subtracting/ adding the term be {T=t . From it, when t&T (and hence e {T=t ?1): which accounts for the hyperbolic shape of the classic Lapicque SD curve.
Originally, eqn. (4) described the SD relationship for extracellular applied current. However, the single-RC equivalent circuit with fixed threshold, where I is the electrode current flowing across the cell membrane can be used with either extra-or intra-cellular stimulation. v~(V {V rest ) is the reduced membrane voltage with V rest the resting value of V : From eqns. (4) and (5), one may also see that b~g(V THR {V rest ), where V THR is the attained membrane voltage at the end of the stimulation (at time T).
Notice that the chronaxie c is not explicitly present in eqn. (4). Notice also that -with very short duration T%t, by the Taylor series decomposition of the exponent (around T~0), one may have either I THR (T)&bt=T or I THR (T)~b½1zt=T: Note that these two different simplifications (and esp. the latter) are 'historical' and depend on which of the two right-hand sides (RHS') of eqn. (4) is used. In the second case only the denominator is developed to first order, while the numerator is truncated at zero-order. The second approximation throws a bridge to Weiss' empirical formula of eqn. (1). I.e. the latter is a simplification of a simplification (i.e. of the 1st-order linear membrane model), capturing best the cases of shortest duration. On the other hand, I THR (T)&bt=T leads to a constant-charge approximation. Interestingly, the latter may fit well also more complex models of the excitable membrane, which take into account ion-channel gating mechanisms, as well as intracellular current flow, which may be the main contributors for deviations from both simple formulas. These 'subtleties' are all clearly described in Lapicque's work, but less clearly by one of the most recent accounts in [16].
Before we continue, it is in order to examine the practical value of numerical optimization to identify energy-efficient waveforms. It is limited for the following reasons. First, it is subject to the rigorous constraints of quantitative equivalence between the model used and the real preparation to which the results should apply. A noteworthy example is provided by the very practice of numerical simulations: often a minute change in parameters precludes the use of a just computed waveform, which is no longer able to elicit an AP in the targeted excitable model. Alas, the same or similar applies hundredfold to the real ES practice.
Second, in the search for minimum-energy waveforms, using numerical mathematical programming algorithms, there is no guarantee about obtaining a globally optimal solution.
Finally, such an approach sheds very little light with respect to the major forces that are at play, and the key factors which determine excitability, such as -for example, the threshold value of membrane potential, whose crossing triggers an AP.
However, the problem at hand is also reminiscent of the search for energy-efficiency in many other physical domains -e.g. ecological car driving. For centuries, physics has tackled similar problems through an approach known as the Least-Action Principle (LAP) [17].
Thus, we first used simple models to derive key analytical results. We then identified generally applicable optimality principles. Finally, we demonstrate how these principles apply also to far more complex and realistic models and their simulations.
The modeling and algorithmic part of this work is laid out in the next section. First, we introduce a simple and general model template. Next we present four most popular specific ionic-current models. Each of these can be plugged in the template to describe an ES target in a single spatial location in excitable-tissue (or alternatively -a space-clamped neural process).
We then examine the conditions for the existence of a finite membrane-voltage threshold for AP initiation. The introduced ionic-current model properties are analyzed to gain important insights into the solution of the main problem at hand.
Two very different ways to identify energy-efficient waveforms are presented in the last two subsections of the Methods. The first relies on a standard numerical optimal-control (OC) approach. The second outlines the LAP in its ES form, which is used to derive a general analytic solution for the energy-optimal trajectories in time of the membrane-potential and stimulation-current.
The Results section presents the model-specific results, applying OC or the LAP. We perform a detailed optimality analysis for both the simple and more realistic models. Comparisons between the two types of approaches, and the quality of their solutions, are made.
Commonly used abbreviations are summarized in Table 1 and symbols -in Table 2.

A General Excitability Model Template
For the equivalent circuit of Fig. 1, I S is the stimulation current. I C is the capacitive current, whose direction is as shown on the Figure when the excitable-membrane's potential is being depolarized. The algebraic sum of all the ionic and all axial currents is represented by I S~IION zI axial , where I axial stands for the algebraic difference (divergence) of in-and out-going axial currents. In the sequel we will use the notation u(t)~I S (t) for the stimulation-current waveform. The latter is our system input, which will be the leverage to refine in order to achieve desirable outcome -reliable triggering of APs in the excitable system. It is customary in the control literature to denote such a signal u(t): Thus, all the currents are linked by the first Kirchhoff circuit law: where -in the most general form, I S depends on membrane voltage V (t) and on the state vector of the ionic channels' gate variables. Unless ambiguous, below we will simplify notation by writing I S (V ): C m (typically around 1 mF =cm 2 , [18]) and V (t) (in mV 's) are the excitable-membrane's capacitance and potential. Equation (6) can be rewritten as: Clearly according to eqn. (7), an outgoing total ionic current opposes the effects of cathodic stimulation, since not all of u(t) is employed toward the main goal of maximizing the V (t) growth, which the reader may have also already deduced from the equivalent circuit of Fig. 1. Conversely, ingoing total current assists the effects of stimulation. Hence, in such a case u(t) may be lower

Specific Single-compartment (Space-clamp) Models
The models here are zero-dimensional (0D). Their spatial extents are confined to a point. This may be contrasted to the multi-compartment cable-like models that we will discuss later, and whose spatial structure is one-dimensional (1D) -i.e. homomorphic to a line.
For single-compartment models there are no axial currents. Hence, I S~IION .
g m is the excitable-membrane's resting (V~V r~2 70 mV ) conductance -in milli-Siemens per unit membrane surface area -e.g. 1 mS=cm 2 . Substituting I ION (V ) from eqn. (8) into eqn. (6) yields a linear first-order model with t~C m =g m~Rm C m the familiar expression for the time constant of such a dynamic model. This model predicts a reasonable resting t& 1 ms. As pointed out in the introduction, this type of model was extensively used by the ES pioneers, [12]. They were particularly concerned with the derivation of analytic expressions for the experimentally observed strength-duration (SD) curves. The latter describe the threshold (minimal) current strength (I THR ), which if maintained constant (i.e. through a rectangular waveform) for a given duration T is likely to elicit an AP in excitable-tissue (see the introductory section).
Even if it may account for a significant part of the sub-threshold variation of the membrane's potential, the linear model lacks a paramount feature -it cannot fire AP's as the latter are due to the highly nonlinear properties of the excitable-membrane's conductance around and beyond the firing threshold.
The Hodgkin-Huxley-type model (HHM). Hodgkin and Huxley (HH) not only proposed a novel way to model ionicchannels but also introduced ionic-channel-specific parameters to fit experimental data [19]. Since, HH-type models have been proposed for many ionic-channels for cardiac to neuroscience applications.
We present one such model from the literature - [20], which has been used to fit experimental data from the central nervous system and particularly the neocortex.
See Tables 3 and 4, which define all the model's variables and parameter values. We consider specifically the Na z v1:6 sodium channel subtype, to which the axon initial segment (AIS) owes its higher excitability [20,21].
The dynamics of a gate-state variable x(t) (where x(t) stands for one of m(t),h(t),n(t)) are described by: Eqns. (6), (9) and (10) define a system of four coupled ODE'swith respect to the four dynamic variables ½V ,m,h,n(t).
Further simplification may reduce the model complexity, maintaining only V (t) as the single dynamic variable. Gatevariable states are factored out by introducing appropriate nondynamic functions of the membrane potential. E.g. in eqn. (9), the fast m gates may be assumed to reach instantaneously m ? (V ), while the far slower h and n gates remain at their resting values (corresponding to a membrane at its resting equilibrium potential V r ).
This model [22] has a second-order nonlinearity, compared to its predecessor -the BVDP model [23], which contains a cubic nonlinearity. The IM will therefore not auto-limit. As in the BVDP, there is a slow second dynamic variable w(t) called the 'recovery current' and its dynamics is described by: The IM responds to supra-threshold stimulation with a wide variety of AP-firing patterns, depending on the particular choices of parameters. Interested in the sub-threshold regimen, we have chosen the ''Spike Latency'' set: b~0:2,c~0:02 [24]. Hence, t w~1 =c is equal to 50 ms. At the time-scale of a single stimulation pulse (lasting at most a few milliseconds), w is virtually a constant.
Here, it may be important to remind the reader that the state of simplest models like the IM needs to be artificially reset after an AP event. However in more complex models (e.g. the HHM), channels that are responsible to revert the system to its resting potential will have a significant effect on the optimal waveform. We will see this in more detail in the results section.

Multi-compartment Models
To expand the scope of our analysis and the applicability of its results, it is essential to also address models of AP initiation and propagation along spatial neural structures. A popular example is the McIntyre, Richardson, and Grill model (MRG 0 02). It was originally used to simulate the effects of ES in the peripheral nervous system and specifically the myelinated axons that form nerve bundles [25]. An adapted version of the same model was recently used to simulate the effects of DBS [7].
Myelinated axon has been pinpointed as the most excitable tissue with extracellular stimulation [26][27][28]. Therefore models like the MRG'02 are of particular interest. Moreover, this model facilitates the illustration of optimality principles as it has only one excitable compartment type -the Ranvier-nodes (RN). The paranodal and other compartments that form the myelinated internodal sections are all modeled as a passive double-cable (due to the myelin sheath that insulates the extracellular periaxonal space) structure, see  Maximum (*2 conductances, in mS=cm 2 : g Na Na z conductance 300 Currents, in mA=cm 2 : Notes: (*1 Typical values are for the Na v1:6 model, [20]; see also Table 4. (*2 These are dependent on (grow with) temperature, the values listed are for T~23 0 C.
(*3 Membrane voltage is either at its resting value V rest ; is depolarized (grows due to stimulation and/or activated sodium Na z ion channels); is repolarized (decays back to V rest , due to the potassium K z ion channels).
(*4 Ionic currents depend on both the membrane voltage and the dynamic state of the ion channels' gates. See Table 4.
The RN compartment is a model of the HH-type: Here two different Na z ion channel subtypes are modeled (please see Table 5 for all the details). The fast subtype (with maximum conductance parameter g Na,f ) is controlled by the opening m and closing h gate states. The persistent subtype (with maximum conductance g Na,p ) is controlled by the p gates. As its name suggests, it has no gate-inactivating states and is non-inactivating. In addition, this model has very slow s gates, associated to its K z ion channel and very fast m gates.
Below we call a fixed point (FP) every V FP value s.t.
The nonlinear dynamics behavior of the RN compartment taken in isolation is quite unlike that of the specific singlecompartment HHM example we provided above. None of its four FPs are stable. Around its unstable 'resting' state (V r = 280 mV ), the zero-dimensional RN's of MRG'02 model yield depolarizing , where T 0~2 3 0 C.
(*2 For a given gate type y of the K z and Na z v1:6 ionic channels, the fractions of open and closed gates are given by the general (Boltzmann-Energy like) template formulae: Thus, the corresponding rates of opening da y =dw and closing db y =dw are sigmoidal functions of w s.t. The actual position of the inflection point (w~0) is determined by the V 1=2 parameter. For the m and n gates, by the l'Hospital-Bernoulli rule, it can be seen that at V m~V1=2 , the opening or closing rates attain half of their max or min, respectively. (*3 For the inactivating gate h of the Na z v1:6 ionic channel: ionic current. I.e. not only does I ION not resist moving away from the resting state, but it actually contributes to automatic firing, with or without any external current! The addition of the passive myelinated spatial structures around the RN's makes the resting state stable, and the problem at hand (of identifying the LAP-optimal ES waveforms) tractable only within a spatial structure. However, this also comes with bonuses. First, the active-passive association brings a very clear-cut picture of the factors at hand that influence AP initiation and propagation. Second, the myelinated double-cable has a very low spatial constant, which provides for a straightforward extension of the single-compartment analysis.
Namely, consider the second term in the more general expression for I S~IION zI axial in eqn. (7). Since around the resting state I ION is always there as a depolarizing factor, it is I axial that needs to be closely considered, see Box in Fig. 2.
The numerical results presented for the MRG'02 in the literature [7,8] often target the mid-cable (center) RN in their ES simulations. This motivated us to use of the method of mirrors to double the model's dimensions at the same computational cost. We consider a long axon (with 41 RN's), which has a relatively low length constant (l 2~1 =(g a r a )). See also Tables 4 and 5. For the RN's l = 167.5 mm vs respectively 2129.7 and 443.2 mm, for the myelinated and the MYSA (paranode) sections. These are paired to significant differences in the passive membrane time We studied extensively all the published accounts of the MRG'02 model and its use for ES modeling [7,8,25]. We also carefully compared parameter values (see Tables 5 and 6) to the ones in the official NEURON models database (senselab.med.yale.edu/modeldb/ShowModel.asp?model = 3810).
Our model implementation originally for [29,30] was done in Matlab (the Mathworks, ver. 7 and above). The code uses CVODES (the Lawrence Livermore National Laboratory, Release 2.7.0) to reliably and robustly solve the related multidimensional system of ODEs. The implementation was validated through extensive comparisons and personal correspondence with the authors of the original model -W.M. Grill [31] and A.G. Richardson, regarding specifically the mismatch between the 2002 publication and its NEURON implementation.

Preliminary Analysis: On the Existence of the AP-firing Threshold
The above ionic-current descriptions differ largely in form and complexity. Yet each of them is capable of capturing some of the essential dynamics properties of excitable living tissues.
In order to elicit an AP through electric stimulation, the membrane's potential V (t) needs to first be driven (depolarized, _ V V w0) to some threshold value V THR , beyond which assisting ionic channels are massively engaged to produce the AP upstroke without the need of any further ES intervention. From eqn. (6) in order to do so, the stimulation waveform needs to be positive and superior to I S (V ,x) at most times -i.e. u(t) needs to overcome the opposing currents. A V THR value is hiding inside each of the above nonlinear flavors of I S (V ,x). Predictably, it is easiest to find the V THR value associated with the IM. Above we saw that the variable w in the IM reacts slowly to changes in V . Hence, one may approximate it by its value at rest: w r~b V r . The resting membrane potential V r is then obtained from the condition I ion,0 (V r )~0, where the subscript 0 indicates that we have assumed w(t)~w r .
The resting potential V r is one of the zeroes of the 2nd-order polynomial in V (t), which characterizes the ionic current. The second zero is V THR . Beyond this threshold the total ionic current Hence, V r = 270 mV and the resting threshold is V THR,0 = 2 55 mV.
We will utilize this simple nonlinear model to complete the picture. If w(V ,t)ww r -i.e. the membrane is not at rest, the point where the total ionic current I ION (V ) switches sign is shifted rightward toward a higher V THR value. For example, for very long durations T??, w?bV : The subscript ? indicates that we have assumed w(t)~bV (t). Predictably, this does not affect the resting potential, since I ion,? (V r )~I ion,0 (V r ). However, V THR,? = 250 mV is higher than the resting threshold V THR,0 .
This reflects the lowering of excitability shortly after an AP, and once the post-AP membrane re-polarization takes place. This is known as refractoriness, which can be either absolute -i.e. no AP can be elicited regardless of how large the stimulation, or relativei.e. larger stimulation current is required -to reach a higher threshold V THR .
Some models of the HH-type have even more complex I ION,? (V ) and thence V THR behavior. This complexity is due to the multiple gate states, which may have very different time constants and hence reach their asymptotic states at different times. In addition, the HH models involve inactivating sodium (Na z ) channels. Hence, excitability may be conditional on attaining the firing threshold within a specific time window. Then V THR may exist only with durations %?. Hence, even over arbitrarily long duration, an arbitrarily low (non-zero) current may never elicit AP's, and may also damage the tissues and the electrodes as irreversible chemical reactions take place.
So, wide stimulation pulses lasting well over some critical duration T CR may not be able to elicit any AP. This is due to the comparable temporal scales of duration T STIM and the time constant t ion of the closing gates associated with depolarizing ionic currents and of the opening gates associated with re-polarizing currents.
Therefore, let us assume that the excitable-membrane's potential is at its resting value V r . Hence, in principle an action potential (AP) can be elicited by stimulation of the fixed duration TvT CR . Therefore stimulation takes place over a finite timehorizon.

Finite-Horizon Optimal-Control (FHOC)
In this approach, the current waveform is the unknown system input signal complying with specific optimality criteria. The optimal pattern u Ã (t) for t[½0,T is sought as a solution of the following constrained minimization problem: where L and R are the constant lower and upper bounds on the values for each u(t) sought. The computational model's dynamical system is introduced in the optimization problem of eqn. (16) in the form of a set of equality constraints. The vector function F(x,u)[R n describes the dynamics of the array of system state-variable trajectories x i (t),i~1 . . . n, resulting from given initial state X(0) and control signal u.
The example developed in the Results section uses the Izhikevich model -eqns. (6) and (11) -with n~2.
The minimized functional, contains the integration term f 0 (X,u) and a final-time (also known as penalty) term H(X(T)) -pulling toward the desired final state X Ã (T). The specific f 0 expression yields minimum electric stimulation power: The penalty term is a convenient way to express the desirable stimulation's outcome -the membrane voltage reaching some predefined threshold-level V THR : Using a general constrained parametric optimal-control approach (e.g. [32]), the objective and equality constraints in eqn. (16) are combined into the Lagrangian: where (t) are the Lagrange multipliers, associated to each of the n equality constraints in eqn. (16) and ( : ) 0 stands for the vector-matrix transpose operator. H~f 0 (x,u)z 0 F(X,u) is known as the Hamiltonian.
The necessary conditions for optimality require that all partial derivatives of the Lagrangian by the system states vanish at the optimal solution to the problem of eqn. (16) -i.e.:
This development is known as mathematical sensitivity analysis and its main purpose is to reveal the impact of a given system parameter (such as u(t) or its initial state X(0)) on the resulting dynamics.
From eqns. (19) and (20): Notice that eqn. (21) describes the adjoint dynamic system iterated in reverse time with a terminal condition provided by the derivative of the h(X(T)) term. To solve the ODE system of eqn.
where Dt is the sampling time, u k~u (kDt) and Moreover, one sees from eqn. (19) that the array (0) is the sensitivity (i.e. the gradient) w.r.t. initial state X(0), i.e.:  A boundary-value problem (BVP), with known initial conditions for X(0) and terminal conditions for l(T), is solved numerically. However, it should also be noted that such solutions may also converge to shallow local minima. For example, the Newton search is guaranteed to produce the 'true' solution when the problem at hand involves a quadratic cost. Here the objective function not only may be non-quadratic, but also may be nonconvex in some manifolds of its high-dimensional parametric space.
Above we described the continuous-time FHOC. The CVODES toolbox readily provides adjoint sensitivity analysis (ASA) capabilities. FHOC is one of the common applications of the latter. Analogously, a discrete-time version may be formulated and solved (see the Results section, where a specific example is developed).

Solving the Problem Analytically: The PLA in ES
Through calculus of variations, here we establish a general form for the energy-optimal current waveform u Ã (t). This approach applies the Principle of Least Action to ES.
Let us assume that T%t ION , where t ION is the time-constant that determines the behavior of the slow gate states of the modeled ionic-channels. Hence, the fast gate states may be approximated by their asymptotic values x ? (V )~lim t?? x(tjV ), while the slow gate states -by their resting values x 0~x? (V r ).
Then an AP can readily be evoked by stimulation from the resting state, and the threshold potential V THR to reach at time T is finite and assumed (without loss of generality) to be known. The energy-efficiency of driving the excitable-tissue membrane potential V (t) from its resting value V r to V THR through a stimulation of fixed duration T satisfies: Since from eqn. (6), u(t)~C m _ V V zI S (V ): As done in the calculus of variations let us perturb the energyoptimal time-course V Ã (t) by the infinitesimal perturbation Eg(t), where g(t) is an arbitrary function of time and E is an infinitesimal scalar.
The necessary condition for S(V ) to have a minimum at E~0 for any g(t) is: To deal with the u Ã (t) _ g g term of eqn. (28), it is integrated by parts : Since the perturbation g(t) respects the boundary-value problem (BVP) with known initial and terminal conditions for V Ã (t) -i.e. g(0)~g(T)~0, then the first RHS term above vanishes. Hence, the only way that eqn. (29) will hold for any g(t) is that we have the Euler-Lagrange-type equation: Equation (30) can also be attained directly using the continuous version of the standard OC formalism [32] (please see also the just presented FHOC subsection above).
Here the Hamiltonian is. The necessary conditions for optimality require that.
Equation (34) is a rather simple system of ordinary differential equations (ODE) that can readily be solved for a given current model I S (V Ã ) to compute the energy-optimal membrane voltage profile V Ã (t). The energy-efficient current waveform u Ã (t) is then computed from eqn. (6).
In the Results section below we illustrate the use of eqn. (34) with several frequently encountered current models.

Results
Here, we first derive some key analytical results using the simplest and clearest models. We then identify generally applicable optimality principles. Finally, we demonstrate how these principles apply also to more complex and realistic models and their simulations.
Part I -Specific Point-model Results, Applying the LAP For the zero-dimensional (single-compartment, space clamp) models introduced in the Methods, here we describe the LAPoptimal waveforms V Ã (t) and u Ã (t), stemming from the general (model-independent) LAP result of eqn. (34).
These simple cases readily illustrate some rather key optimality principles resulting from a LAP perspective. We will discuss these optimality principles as we go, and will summarize them at the end of this subsection.
Linear sub-threshold model. Replacing I S (V Ã ) in eqn. (34) with I ION (V ) from eqn. (8): t~C m =g m~Rm C m is the membrane's time constant and for expediency V :V Ã and V r = 0. The general solution of eqn. (35) is: Given the boundary conditions V (0)~0 and V (T)~V THR : A result similar to eqn. (37) is obtained by [33], using a slightly different (less direct or general) optimal-control approach.
From eqn. (37) one can see that V Ã (t)=V THRs inh (t=t)= sinh (T=t)~t=T -i.e. it has a linear rise, especially with T%t. Here T = 100 ms and t = 1 ms (computed using typical values from the literature for g m = 1 mS=cm 2 and C m = 1 mF =cm 2 ). Figure 4 presents the LAP energy-optimal stimulation profiles V Ã (t) and u Ã (t) for a short and a long stimulus duration T STIM and three membrane time constant t values.
Before we go on, it is useful to investigate the conditions for a growing exponent (GE) waveform to outperform the SQR waveform.
First, u Ã GE (t) has a very rapid rise. Hence, its optimal duration T Ã GE will be short. Second, it is noteworthy that in [33] t = 30.4 micro-seconds! Hence, injected current rapidly leaks out. However even with the above extreme t value, at its optimal duration T Ã SQR the SQR wave does just 22% worse, which means that the SQR is among the best candidates for its robustly good performance.
Second, in multiple cases, the energy-optimal LAP waveform u Ã (t) looks a lot like a 'classical' rectangular waveform. From eqn. (8), we may also see that, with V r = 0, V THR = 1, the max. value of I ION (V ) is equal to 1 and is attained as the membrane potential reaches the threshold V (T)~V THR . If we then replace  6), we see that a waveform u(t) -that brings V (t) from V r to V THR at a constant rate, is the time-constant waveform u(t)~k Ã~( V THR {V r )=T. For this example, k Ã1 0&I ION (V ), which explains why u Ã (t) is that close to a rectangular waveform.
As a matter of fact, for very short stimulation times, the k Ã tend to be high, while I ION (V ) tends to be linear. Hence, the 'classic' rectangular (or square, SQR) waveform tends to also be close to energy-optimal.
Such facts are rather important as they lead us below (as evidence is accumulated) to a general form not only of V Ã (t), but also of u Ã (t).
Comparative properties the V(t) growth profiles. The GE waveform may be an SQR waveform in disguise. I.e. some linear growth of the membrane voltage may still fit the one obtained upon ES with a GE. The motivation for this is in eqn. (36), where the first term vanishes with T&t.
Finally, the total electric charge conveyed by the ES source may have to be considered. For example, in the LM of eqn. (8) the total charge consists of a capacitive charge to raise the membrane voltage by a given amount (to V THR ), and resistive charge Ð TSTIM 0 V (t)=Rdt. A similar situation occurs in the MRG 0 02 model due to the opposing axial currents.
So let us solve the following auxiliary problem: Find a linear fitV V (t)~max½a(t{b),0 to the growing exponent V (t)~(e t=t {1)=(e TSTIM =t {1), so that the ES source conveys the same resistive charge in the time interval t[½0,T STIM . I.e. we want that: Here, for simplicity (and without any loss of generality) we have assumed V r~0 and V THR~1 .
The latter result promotes intuition: with large opposing currents optimal ES cannot afford to last long. The transition of the membrane voltage from its rest to a threshold value is best performed rapidly. Hence, the shape of the V (t) growth profile depend on the T STIM =t ratio. As seen, for T STIM %t, the optimal u(t) is close to rectangular, while with T STIM &t, the GE is in effect equivalent to doing nothing for at least half of the duration, and then to a SQR waveform of at least doubled amplitude.
With quite similar reasoning, one can demonstrate that a 1storder membrane voltage growth profile V (t)~(1{e t=t )= (1{e TSTIM =t ) in the time interval t[½0,T STIM is suboptimal and equivalent to linear growth, which has about twice longer duration.
Izhikevich model. Replacing I S (V Ã ) in eqn. (34) with the I ION (V ) approximations from eqn. (14) or (15), see Box in Fig. 5: As in the preceding model V :V Ã . Note that the dynamics of eqn. (38) has all FP's of I ION (V ), as well as a third FP at V~0:5(V r zV THR ), contributed by the derivative term I 0 ION (V). Equation (38) can be solved analytically. However, it provides the solution in an implicit form and involves an incomplete elliptic integral of the first kind. Hence, we used the Matlab bvp4c BVP solver with boundary conditions V Ã (0)~V r and V Ã (T STIM )Ṽ THR . Figure 5 illustrates the energy-optimal LAP solution u Ã (t) and the corresponding membrane voltage profile V Ã (t). The I ION,0 (V ) approximation of the ionic current is used for a case of very short duration (T STIM = 10 ms) and the I ION,? (V ) approximation is used for a case of long duration (T STIM = 5 ms).
It is important to notice that -as with the LM model above, u Ã (t)&k Ã zI ION (V ), where k Ã~( V THR {V r )=T (see the Box in Fig. 5).
According to eqns. (14) and (15) the opposing current in the IM can be presented in the general form: To see how the optimal ES is affected by the level of opposing current, it is more than tempting to experiment with different gain values.
Hence, 3 gain cases are plotted in Fig. 5 -for the nominal gain (cyan traces) and two additional cases: the opposing current I S (V Ã ) is either doubled (gain = 2, red traces) or decreased twofold (gain = 1/2, black traces). As could be intuitively expected from the general equation (24), when I ION (V )?0 (very low ionic currents): By the Cauchy-Schwartz inequality in the space of continuous real functions, it is straightforward to show that the voltage trajectory V Ã (t) that minimizes eqn. (40) is such that _ V V Ã (t)~k Ã , where k Ã is determined from the boundary conditions satisfied by V Ã (t). Hence: Just as in the preceding model, it is also V Ã (t)=V THR~t =T STIM with the shorter durations -which justifies the use of the resting approximation I ION,0 (V ).
HHM. Here the I S (V Ã ) of eqn. (34) is replaced with the resting-state -I ION,0 (V ), or asymptotic-state -I ION,? (V ) ionic current approximations (see the Box in Fig. 6).
Toward I ION,0 (V ) the gate-state variables are factored out as follows: The fast state m(t)&m ? (V ), while the slower variables h(t)&h r~h? (V r ), and n(t)&n r~n? (V r ) are approximately at rest, assuming very short durations. Conversely, and assuming very long durations, toward I ION,? (V ) all gate variables are approximately at their asymptotic value, corresponding to a given membrane voltage V (t) (see Methods).
As with the IM, we used bvp4c to numerically solve the BVP of eqn. (34) with boundary conditions V Ã (0)~V r and V Ã (T STIM )Ṽ THR . Figure 6 follows a very similar format to Fig. 5. Similarly to eqn. (39) above, I S (V ) can also be assumed higher or lower. All the maximal ionic conductances in the HHM (see also Table 3) are temperature-dependent and are linearly proportional to the coefficient k T : where Q 10~2 :3 and T 0 = 23uC. Hence with T = 37uC, according to eqn. (42) k T = 3.2094. Let this be our standard case (gain = 1).
As we did with the IM, 3 gain cases are plotted in Fig. 6 for I S (V Ã )~gain|Î I ION (V ). For the two additional cases the opposing current I S (V Ã ) is either doubled (gain = 2, red traces) or halved (gain = 1/2, black traces).
Once again -as with the LM and IM models above, u Ã (t)&k Ã zI ION (V ) (see the Box in Fig. 6).
Numerical model simulation and optimal control. The IM was also evoked in the FHOC Methods section. It is therefore interesting to contrast the results of the LAP and FHOC approaches in identifying energy-optimal ES waveforms for the same ionic current model. For such comparison, the IM has the clear advantage of hiding no implementation specifics inside a black box.
The FHOC formalism (see Methods) is computationally efficient, but it is also subject to the similar limitations as most of the ad-hoc search approaches. Iterative numerical optimization requires an initial guess for the solution, and trying different starting arrays u (0) may alleviate a bit the propensity to converge to shallow local energy-minima.
Here it is also important to realize that in eqn. (16) the two terms to minimize in the F (u) functional (a function of functions), namely the energy cost (17) and the penalty (18) may conflict each other. When the penalty gain K penalty in (18) is too low, the search will identify a lower-energy solution u, which however does not bring the membrane potential V k up to the desired threshold value -i.e. V M %V THR . Conversely, a too high penalty gain K penalty will identify a very high-energy solution u, which is not only costly, but the membrane potential may also overshoot the threshold, since The linear-ramp voltage profile yields the best P performance for most of the durations. As in Fig. 8 notice that the P and Q values are quite similar for the linear and exponential cases, for T STIM respectively 2 and 5 ms; and also for the 1 st -order and linear cases, for T STIM respectively 0.2 and 0.5 ms. Toward the P values electrode impedance of 1 MV is assumed. Contrasted: SQR stands for the square (or rectangular) stimulation waveform. doi:10.1371/journal.pone.0090480.g010 the 'getting there' is underestimated for the sake of the very last simulation steps.
As seen from Fig. 7 Panel B (which uses the I ION,? (V ) approximation of the ionic current for the relatively long duration T STIM = 2 ms), the linear growth profile is a reasonable estimate for the optimal membrane voltage profile V Ã (t). Hence: where k Ã is given by eqn. (41). When u (0) is close to the LAP estimate u Ã (t) of eqn. (43), the FHOC iteration also consistently ends close to there (see Fig. 7, panel B). The cyan traces on Fig. 7 are the u Ã (t) and the resulting V Ã (t). With the LAP estimate, the FHOC approach resulted in a final membrane potential reasonably close to the desired threshold value -i.e. V (T STIM )~{50:106&V THR~{ 50, even if the IM was simulated with the discretized LAP waveform u Ã (t) (Dt = 10 ms).
The black traces illustrate the FHOC solution, computed for two different u (0) choices. For Panel A, u (0) was chosen to be all zeros. When all time-step entries u (0) were chosen to be equal to the upper bound U = 30 (data not shown), due to the (discontinuous) AP event occurring mid-way the temporal horizon, the Matlab's fmincon solver remains stuck to the initially provided values.
Except for the case in Panel B, the K penalty meta-parameter had to be kept high (K penalty = 70) in order to respect the terminal constraint of V (T STIM )&V THR .
The total energy costs (all expressed as 2-norms of the obtained best u(t)) are respectively 161, 153.2 and 423.4 (for the discretetime version) 186.7, 159.1 and 334.2 (for the continuous-time version).
Comparing these to P(u Ã ) = 153.2 (discrete-time) and = 157.4 (continuous-time), the LAP-based solution is comparable to or superior than the FHOC solutions. The numerical FHOC solution on Fig. 7, panel A has converged to a local extremum. Note that a post-hoc correction (simple DC offset) is applied to the LAP-based estimate, which adjusts for the overshoot of V THR when simulating the full (two-dimensional) IM. The overshoot is due to using the one-dimensional approximation, eqn. (15).
The results obtained here nicely illustrate multiple aspects of identifying energy-efficient waveforms through numerical model simulation and optimization. Clearly, pairing theoretical insights with numerical tools carries the best success potential.

Part I Results Summary
A number of more general observations on u Ã (t) can be made looking at the results this far.
Probably, the most significant result is that the use of LAP reduces the problem to the BVP, defined by eqn. (34), with V Ã (0)~V r and V Ã (T STIM )~V THR . We still need to have a very good idea of both I S (V ) and V THR to successfully solve for V Ã (t), and thence for u Ã (t), in a given particular situation.
We identify also the following key and practice-oriented optimality principles resulting from the LAP perspective.
1. The optimal sub-threshold membrane potential growth profile with relatively short durations T STIM and low membrane conductivity: First, in all simple models we used up to here, the solution V Ã (t) of the ODE system, defined by eqn. (34), is quite close to a linear growth from V Ã (0)~V r to V Ã (T STIM )~V THR . Second, with the total current I S (V )&0 (e.g. low leak), then from eqn. (6), it follows that u(t) will be exactly proportional to the rate of change of the membrane's potential V(t). If _ V V Ã (t)&const, then u Ã (t) is close to a SQR waveform.
2. The energy-efficient waveform depends directly on the temporal shape of currents at the AP initiation site. 3. The targeted V THR membrane voltage threshold depends on stimulation duration, with a tendency to increase with T STIM . 4. The exponential growth membrane voltage profiles V (t) are equivalent to linear growths of shorter duration.

Part II -Multiple-compartment Model Results
Here we first extend the general (model-independent) LAP result of eqn. (34) to spatial-structure models (non-zero-dimensional, multi-compartment), which involve membrane-voltage distribution and propagation along cable structures.
LAP result generalization to multi-compartment models. There is a combinatorial explosion in both the number of parameters and the number of ways that multi-compartment models can be put together and used. Hence, there is much more than one way of generalizing the LAP result of eqn. (34).
Here we briefly present a variant, which appears to be one of the most straightforward generalizations. With a multi-compartment model, eqn. (7) can be rewritten as: Without loss of generality, we used the variable Z to represent any 'spatial' model dimension. It could even stand for the compartment index in a discretized implementation. Now, eqn. (7) is a partial DE, depending both on the temporal and the spatial model dimensions.
Assuming that we are free to manipulate u(t,Z) in every compartment as we wish, the derivation sequence from eqn. (23) to eqn. (30) (see the LAP subsection in the Methods) still applies yielding a family of equations 'parameterized' by the location coordinate Z.
Hence, we may obtain the generalization of eqn. (34) as: Like the extended eqn. (44), eqn. (45) is a partial DE, depending on both temporal and spatial boundary conditions. In particular, V THR becomes a function of Z. It is no longer a single variable, but a whole spatial profile, subject to conditions such as the safety factor for propagation introduced in the cardiac literature [34].
The MRG'02 model: Toward upper bounds on V THR (T STIM ). Multi-compartment models add complexity unseen with the single-compartment models. Wongsarnpigoon & Grill [8] used the peripheral-axon MRG'02 model [25] in a genetic-programming search for energy-efficient stimulation waveforms. The approach was somewhat similar to the FHOC described above. After thousands of iterations simulating the MRG'02 model, the identified waveforms were reminiscent of noisy truncated and vertically offset Gaussian's ( Fig. 2 in [8]). In the light of analysis this far one might think that this reflects the shape of I S (V ) for V ranging from the resting value (280 mV) to some threshold V THR .
In this work stimulation is assumed to be intracellular and at just one spatial location (Z~0, the center RN, see Methods) along the cable structure.
To suggest a version of optimal waveforms u Ã (t) for the MRG'02 model, we first estimate the membrane voltage threshold for each duration. One analytic way toward such estimates is  Energy-Optimal Electrical-Stimulation readily provided by the MRG'02 model. Recall also that with simpler models V THR showed a tendency to increase with T STIM . Figure 8 presents a family of ionic current I ION (V ,Z) approximations at the target site (Z~0), for a set of durations T STIM . For each of the durations we assume that the membrane voltage trajectory V (t) evolves according to a linear ramp from rest V r to threshold V THR . As the latter is unknown, we produced one such ramp for each V value on the horizontal (independentvariable) axis of the figure, and then computed the corresponding ionic current I ION (V ) as described next.
Toward gross estimates of V THR , we first solve approximately eqn. (10) for each gate-state: where x 0 is the gate-state value at rest and V V~(V r zV )=2 is the average excursion from the resting membrane voltage. Figure 8 shows the obtained approximate ionic currents I ION (V (t)) as a function of just V for three very different durations -T STIM = 0.02, 0.5 and 5 ms. For T STIM = 5 ms, the Box in the same figure illustrates the estimated proportions-to-rest x x(T STIM jV )=x ? (V r ) for each of the 4 gate-state variables, at the end of stimulation.
Why does such an analysis provide upper bounds on V THR (T STIM )?
First, from the Box of Fig. 8 we can see that indeed the dynamics of the fast Na z ion channel subtype evolves before that of the other ion channels. Particularly, we see that the estimate for inactivating h gates suggests they are completely closed for T STIM = 5 ms and once V reaches around 240 mV .
On the other hand from the main Fig. 8, one can see that this analysis gives the intervals V[½V rest ,V UB in which the approximate ionic currents I ION (V )v0 (i.e. remain depolarizing).
Clearly if V THR (T STIM ) is not reasonably within ½V rest ,V UB , no miracle would yield an AP at the target location, since I ION becomes repolarizing outside of these bounds.
Interestingly, the analysis also predicts lowering of V THR with longer durations. This result is exactly the opposite of what was observed with the simpler models of the HH-type, where I ION was repolarizing for V [½V rest ,V THR .
The numerical experiments we conducted were fully consistent with the above predictions, and some upper bounds were also quite tight.
The MRG'02 model: numerical experiments. We conducted four series of numerical experiments in search of the optimal waveforms u Ã (t) for the MRG'02 model. Each series was computed for the same set of 9 durations T STIM = 20, 50, 100, 200, 400 and 500 ms; 1, 2 and 5 ms (for the sake of better visibility, only the most representative subsets are illustrated in full detail).
The four series differed by the chosen voltage-clamp temporal growth profile V (t,0) at the targeted RN location and A baseline series involved finding the threshold rectangular stimulation amplitude. In all series, the constraint was to observe a propagating AP at the latest within 1 ms after the end of stimulation.
With DV~V THR (T STIM ){V r , where the minimum V THR was found (with 0.001 mV tolerance) using the same type of goldensection search algorithm as per the optimal SQR amplitude.  Figure 13. Propagating AP due to an optimal SQR (rectangular) waveform, T STIM = 100 ms: For the shortest durations, the plain rectangular waveform outperforms by P the ones associated to the linear-ramp voltage profile. One can see clearly that the steep rise of the SQR waveform yields an early superlinear ramping of the membrane voltage. However, the rectangular waveform requires a lot more charge Q to be transferred (see Fig. 10). doi:10.1371/journal.pone.0090480.g013 And the three LAP-driven series were: linear growth exponential growth 1-st order growth The corresponding u(t,0) ES waveforms were computed from eqn. (44) with Z~0.
The MRG'02 model: numerical results. Figure 9 and table 7 illustrate the obtained V THR as a function of T STIM .
The computed optimal values of V THR are often similar for two adjacent durations either between the linear and 1-st order, or between the linear and exponential growth (EG). 1-st order is usually similar to its right-hand linear neighbor (for the next longer duration). Conversely, EG is similar to its left-hand linear neighbor (for the previous shorter duration). This is consistent with and best interpreted in the light of our growth-profiles comparison (see the dedicated subsection on page 13). There we saw that indeed an EG V (t) trajectory is approximately equivalent to linear growth of about twice shorter duration. As for 1-st order growth, clamping the voltage to its plateau will tend to be similar to a linear growth of about twice longer duration. Recall also that 1-st order is the 'reverse-time' analog of EG. Figure 10 and tables 7, 8 illustrate the obtained optimalwaveforms' energy P and charge-transfer Q values as a function of T STIM .
The linear-growth strategy is the one that tends to perform best across the board, except for the 2 longest durations, and as predicted by the comparative (linear vs exponential growth) analysis, based on the 0D LM. Figure 3 illustrates the propagating AP's, corresponding to the two representative linear and exponential voltage-clamp temporal growth profiles at the stimulation site V (t,0). The figure also shows the spatial profiles of the membrane voltage and intracellular potential at the end of stimulation for the two growth cases.
Consistently with the analysis in the subsection on the comparative properties of the V (t) growth profiles, we found out that the spatial distributions of membrane voltage and intracellular potentials at the end of stimulation were reasonably similar -e.g. between the optimal linear growth voltage-clamp for T STIM = 2 ms, Fig. 3 (Panels A, C) and the optimal exponential growth with T STIM = 5 ms, Fig. 3 (Panels B, D).
Note that we expect from an approximately globally optimal stimulation waveform u Ã (t) to yield a specific distribution of membrane voltages V (T STIM ,Z) at the end of the stimulation. We call this distribution tentatively the invariant spatial profile of the membrane voltage. Importantly, such a profile will differ for any different duration T STIM even when the corresponding waveform u Ã (t) is globally optimal. This is due for example to the small spatial constant l, which controls the spatial diffusion with time.
However, if the spatial profile is about the same for different durations T STIM and the corresponding different waveforms u Ã (t) (see Panels B and D in Fig. 3), then both waveforms may be optimal. Recall that linear fits to both the optimal 1-st order growth and the optimal exponential growth with durations T STIM = 5 ms have duration &0:46|T STIM = 2.3 ms. Thus, all of the above cases may yield quasi-invariant spatial potentials at the end of stimulation, and may also be otherwise similar.
For two representative linear-growth cases Fig. 11 illustrates the corresponding waveforms u Ã (t) and their construction in detail.
Finally, Fig. 12 uses the same-vertical-scale to compare the relative contributions of the growth rate and the compensated repolarizing node currents for each different duration. The waveforms' offsets (due to k Ã ) are inversely proportional to duration. This readily compares qualitatively with the results in [8]. Especially for very short durations (e.g. T STIM~2 0ms), the optimal waveform u Ã (t) has a significant rectangular component (see also the optimality-analysis for the simple 0D models). Further parallels may be made for the relatively shorter durations (ƒ200ms).
Numerous essential differences in the approach preclude further objective comparisons. Interestingly however, for the longer durations (T § 0.5 ms) the results in [8] show very little (if any) variation with T STIM (there called pulse-width, PW).
Finally, with long PW's in [8] most of the stimulation's energy is delivered toward the middle of the active period. This late and peaky delivery requires additional analysis and comparisons of the actually achieved waveform-energy levels, which cannot be done in its details at this time. However, we return to the late delivery policy in the Discussion (see below), where it is deemed equivalent to a shorter-duration case.
The latter provides a clue why such significant delivery differences would not be at odds with the very narrow 95% confidence intervals that resulted from the genetic algorithm in [8], and seeming to preclude different optimal waveforms.

Discussion and Conclusions
In eqn. (23), we addressed directly the electric power required for driving the excitable-tissue membrane potential V (t) from its resting (V r ) to its threshold value (V THR ) through a stimulation of fixed duration. Through the LAP perspective, we obtained eqn. (34) -a general (model-independent) description of the energyoptimal time-course of the excitable-tissue's membrane potential V Ã (t).
We would like to bring the reader's attention to three specific conclusions.
The first is related to the intuition gained with respect to the evolution of the membrane potential V Ã (t). This optimality principle is best demonstrated by the simplest linear sub-threshold model (LM). Let ES circumstances be characterized by large opposing currents (e.g. the leak LM current) over long durations. This situation is physically analogous to filling with water a bucket which has large holes in its bottom. Since only the final outcome is important (i.e. we want the bucket full at the final time T), the best policy is to do nothing for most of the duration and then be able to dump a very large amount of water in the bucket over very short time. From experience, we know that works for even an unplugged sink. Moreover, we saw that the same intuition transfers to more refined models (e.g. the HHM or the MRG'02) as do nothing for most of the duration means that we are still around the resting V and hence there is no danger of Na z ionic-channel deactivation.
The second take-home message is that the use of LAP principles jointly with numerical approaches (e.g. the classical FHOC) provides a mathematically sound and practical waveform optimi-zation approach, providing more assurance toward the quality of the final outcome.
And finally, a note of humility is in perfect order. In this work we just slightly opened the door to using the LAP ideas for optimal ES. There are many more aspects to tackle than the ones that we can address in this short paper as 'proof of concept'. In particular we would like to extend the method for extracellular stimulation in forthcoming work. The motivation for doing so is at least twofold. On the one hand, extracellular stimulation has far more practical relevance. On the other hand, the only way we could rigorously employ the general LAP solution of eqn. (45) is to consider a model where we are free to manipulate u(t,Z) in every compartment or at every spatial location.
A direction for such manipulation is provided by the activating function concept [15,20,25], which supplies every compartment with a virtual injected current. In the context of extracellular stimulation, we will also have to properly address the conditions for stable AP propagation (see [15,35] for an extensive treatment of the subject). The optimal pattern of extracellular potentials (size of depolarized and hyperpolarized regions) depends on the distance to the electrode. These conditions would also naturally provide the spatial voltage profile at the end of the stimulation, needed to properly solve the PDE of eqn. (45).
Here we took a shortcut path by assuming that intuitions gained with single-compartment models suffice. This may be partially true with the specific MRG'02 setup that we addressed, but does not hold in general. Hence, the LAP results are approximate. A clue is provided by the slightly lower P values of the optimal rectangular waveform, for T STIM = 100 and 200 ms -see Table 9. As can be seen from Fig. 9, no benefit in terms of lower V THR can be associated to the steep rise of the rectangular waveform, since V THR is expected to be higher, esp. for dramatically shorter durations. This was further confirmed by numerical testing with dual linear (high/low rate) V(t) rise schedules (data not shown), which all had inferior performance to the baseline simple lineargrowth protocol. However, the rectangular waveform also leads to steep capacitive decay of V (t) at the end of the stimulation, which may trigger specific patterns of additional depolarizing currents.
For the shortest durations, the plain rectangular waveform outperforms by P the ones associated to the linear-ramp voltage profile (see Fig. 10). On Fig. 13 one can see that the steep rise of the SQR waveform yields an early super-linear ramping of the membrane voltage. However, the rectangular waveform requires a lot more charge Q to be transferred.
In practical situations many more additional aspects need to be addressed. E.g. stimulation needs to be charge balanced. This is a necessity for implanted devices and also debatably important for transcutaneous applications. Such stimulation will have an effect on the optimal threshold intensity of the cathodic pulse [36]. One would expect that a pre-or post-anodic pulse would also have a significant effect on the optimal waveform. Moreover, its own shape would be subject to optimization -e.g. to minimize the overall energy level required -a cost suitable for the design of implanted devices.
We hope that the analysis and numerical evidence provided in this work may convince the reader of the practical benefits of applying the LAP principles toward the design of energy-efficient ES.