Corrections
10 Apr 2015: The PLOS ONE Staff (2015) Correction: Energy-Optimal Electrical-Stimulation Pulses Shaped by the Least-Action Principle. PLOS ONE 10(4): e0124704. https://doi.org/10.1371/journal.pone.0124704 View correction
12 Dec 2014: The PLOS ONE Staff (2014) Correction: Energy-Optimal Electrical-Stimulation Pulses Shaped by the Least-Action Principle. PLOS ONE 9(12): e115928. https://doi.org/10.1371/journal.pone.0115928 View correction
Figures
Abstract
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.
Citation: Krouchev NI, Danner SM, Vinet A, Rattay F, Sawan M (2014) Energy-Optimal Electrical-Stimulation Pulses Shaped by the Least-Action Principle. PLoS ONE 9(3): e90480. https://doi.org/10.1371/journal.pone.0090480
Editor: Dante R. Chialvo, National Research & Technology Council, Argentina
Received: November 11, 2013; Accepted: January 30, 2014; Published: March 13, 2014
Copyright: © 2014 Krouchev et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: The Fonds de recherche du Quebec - Nature et technologies and the Natural Sciences and Engineering Research Council of Canada provided funding for this work. S.D. was supported by the Vienna Science and Technology Fund (WWTF), Proj.Nr. LS11-057. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
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]–[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]–[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 model-specific. 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]–[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](1)where
is the stimulus duration,
is called the rheobase (or rheobasic current level) and
is the chronaxie.
The most expedite way of introducing the rheobase and chronaxie would be to point to eqn. (1) and notice that:(2)and
(3)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, Both Lapicque’s own writings - [11]–[13], and more recent work are at odds with the linear-charge approximation. Already in 1907 Lapicque was using a linear first-order approximation of the cell membrane, modeled as a single-RC equivalent circuit with fixed threshold:
(4)with time constant
and
are the membrane capacity and conductance respectively.
The second form of eqn. (4) is easily obtained by subtracting/adding the term . From it, when
(and hence
):
which accounts for the hyperbolic shape of the classic Lapicque SD curve.
Originally, eqn. (4) described the SD relationship for extra-cellular applied current. However, the single-RC equivalent circuit with fixed threshold, where is the electrode current flowing across the cell membrane
(5)can be used with either extra- or intra-cellular stimulation.
is the reduced membrane voltage with
the resting value of
From eqns. (4) and (5), one may also see that
where
is the attained membrane voltage at the end of the stimulation (at time
).
Notice that the chronaxie is not explicitly present in eqn. (4). Notice also that - with very short duration
by the Taylor series decomposition of the exponent (around
), one may have either
or
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,
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.
Methods
A General Excitability Model Template
For the equivalent circuit of Fig. 1, is the stimulation current.
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
where
stands for the algebraic difference (divergence) of in- and out-going axial currents. In the sequel we will use the notation
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
is the intra-cellular stimulation current.
is the capacitive current. The direction of the latter is for a case of depolarizing the membrane’s voltage (i.e. the inside of the cell wall becoming more positive). The algebraic sum of all the ionic and all axial currents is represented by
where
stands for the algebraic difference (divergence) of in- and out-going axial currents.
Thus, all the currents are linked by the first Kirchhoff circuit law:(6)where - in the most general form,
depends on membrane voltage
and on the state vector of the ionic channels’ gate variables. Unless ambiguous, below we will simplify notation by writing
(typically around 1
, [18]) and
(in
’s) are the excitable-membrane’s capacitance and potential. Equation (6) can be rewritten as:
(7)
Clearly according to eqn. (7), an outgoing total ionic current opposes the effects of cathodic stimulation, since not all of is employed toward the main goal of maximizing the
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
may be lower than when it is estimated assuming the absence of membrane conductivity. Let us elucidate right away by providing typical examples.
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. homo-morphic to a line.
For single-compartment models there are no axial currents. Hence, .
Linear Sub-threshold model (LM).
(8)
is the excitable-membrane’s resting (
−70
) conductance - in milli-Siemens per unit membrane surface area - e.g. 1
. Substituting
from eqn. (8) into eqn. (6) yields a linear first-order model with
the familiar expression for the time constant of such a dynamic model. This model predicts a reasonable resting
1
.
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 (), which if maintained constant (i.e. through a rectangular waveform) for a given duration
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 ionic-channels 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 sodium channel subtype, to which the axon initial segment (AIS) owes its higher excitability [20], [21].
The dynamics of a gate-state variable (where
stands for one of
) are described by:
(10)
Eqns. (6), (9) and (10) define a system of four coupled ODE’s - with respect to the four dynamic variables .
Further simplification may reduce the model complexity, maintaining only as the single dynamic variable. Gate-variable states are factored out by introducing appropriate non-dynamic functions of the membrane potential. E.g. in eqn. (9), the fast
gates may be assumed to reach instantaneously
, while the far slower
and
gates remain at their resting values (corresponding to a membrane at its resting equilibrium potential
).
The Izhikevich model (IM).
(11)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
called the ‘recovery current’ and its dynamics is described by:
(12)
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: [24]. Hence,
is equal to 50
. At the time-scale of a single stimulation pulse (lasting at most a few milliseconds),
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 (). 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]–[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 Fig. 2.
The RN compartment is a model of the HH-type:(13)
Here two different ion channel subtypes are modeled (please see Table 5 for all the details). The fast subtype (with maximum conductance parameter
) is controlled by the opening
and closing
gate states. The persistent subtype (with maximum conductance
) is controlled by the
gates. As its name suggests, it has no gate-inactivating states and is non-inactivating. In addition, this model has very slow
gates, associated to its
ion channel and very fast
gates.
Below we call a fixed point (FP) every value s.t.
. From eqn. (7) with
,
The nonlinear dynamics behavior of the RN compartment taken in isolation is quite unlike that of the specific single-compartment HHM example we provided above. None of its four FPs are stable. Around its unstable ‘resting’ state ( = −80
), the zero-dimensional RN’s of MRG’02 model yield depolarizing ionic current. I.e. not only does
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 in eqn. (7). Since around the resting state
is always there as a depolarizing factor, it is
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 (). See also Tables 4 and 5. For the RN’s
= 167.5
vs respectively 2129.7 and 443.2
, for the myelinated and the MYSA (paranode) sections. These are paired to significant differences in the passive membrane time constant (
). For the RN’s
= 0.29
vs respectively 20 and 2
, for the myelinated and paranode sections. The cable end-conditions are formed by virtual compartments with membrane at rest
= −80
. This choice is further motivated by the results of model simulations - namely the relatively little spread of potentials at the end of stimulation lasting up to a few milliseconds (see Fig. 3).
The peaks of at the Ranvier nodes are due to the direct exposure to the extracellular medium, which is unlike that of the myelinated sections in the double-cable MRG’02 model.
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 multi-dimensional 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 needs to first be driven (depolarized,
) to some threshold value
, 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
at most times - i.e.
needs to overcome the opposing currents.
A value is hiding inside each of the above nonlinear flavors of
. Predictably, it is easiest to find the
value associated with the IM. Above we saw that the variable
in the IM reacts slowly to changes in
. Hence, one may approximate it by its value at rest:
. The resting membrane potential
is then obtained from the condition
, where the subscript
indicates that we have assumed
.
The resting potential is one of the zeroes of the 2nd-order polynomial in
, which characterizes the ionic current. The second zero is
. Beyond this threshold the total ionic current switches its sign. So eqn. (11) becomes:
(14)
Hence, = −70 mV and the resting threshold is
= −55 mV.
We will utilize this simple nonlinear model to complete the picture. If - i.e. the membrane is not at rest, the point where the total ionic current
switches sign is shifted rightward toward a higher
value. For example, for very long durations
,
:
(15)
The subscript indicates that we have assumed
. Predictably, this does not affect the resting potential, since
. However,
= −50 mV is higher than the resting threshold
.
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 relative - i.e. larger stimulation current is required - to reach a higher threshold .
Some models of the HH-type have even more complex and thence
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 (
) channels. Hence, excitability may be conditional on attaining the firing threshold within a specific time window. Then
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 may not be able to elicit any AP. This is due to the comparable temporal scales of duration
and the time constant
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 . Hence, in principle an action potential (AP) can be elicited by stimulation of the fixed duration
. Therefore stimulation takes place over a finite time-horizon.
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 for
is sought as a solution of the following constrained minimization problem:
(16)
where
and
are the constant lower and upper bounds on the values for each
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 describes the dynamics of the array of system state-variable trajectories
, resulting from given initial state
and control signal
.
The example developed in the Results section uses the Izhikevich model - eqns. (6) and (11) - with .
The minimized functional, contains the integration term and a final-time (also known as penalty) term
- pulling toward the desired final state
. The specific
expression yields minimum electric stimulation power:
(17)
The penalty term is a convenient way to express the desirable stimulation’s outcome - the membrane voltage reaching some pre-defined threshold-level :
(18)
Using a general constrained parametric optimal-control approach (e.g. [32]), the objective and equality constraints in eqn. (16) are combined into the Lagrangian:(19)where
are the Lagrange multipliers, associated to each of the
equality constraints in eqn. (16) and
stands for the vector-matrix transpose operator.
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.:(20)
Here the ‘vector-matrix’ notations or
, where
, mean respectively
or
,
.
This development is known as mathematical sensitivity analysis and its main purpose is to reveal the impact of a given system parameter (such as or its initial state
) on the resulting dynamics.
From eqns. (19) and (20):(21)
Notice that eqn. (21) describes the adjoint dynamic system iterated in reverse time with a terminal condition provided by the derivative of the term. To solve the ODE system of eqn. (21), the achieved forward dynamics of eqn. (16) needs to be already computed.
Similarly, all partial derivatives of the Lagrangian by vanish at the optimal solution to the problem of eqn. (16) - i.e.
:
(22)where
is the sampling time,
and
Hence, eqn. (22) yields all components of the gradient w.r.t. , which enables the use of gradient-based quasi-Newton search routines (e.g. fmincon from the Matlab optimization toolbox).
Moreover, one sees from eqn. (19) that the array is the sensitivity (i.e. the gradient) w.r.t. initial state
, i.e.:
A boundary-value problem (BVP), with known initial conditions for and terminal conditions for
, 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 non-convex 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 . This approach applies the Principle of Least Action to ES.
Let us assume that , where
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
, while the slow gate states - by their resting values
.
Then an AP can readily be evoked by stimulation from the resting state, and the threshold potential to reach at time
is finite and assumed (without loss of generality) to be known. The energy-efficiency of driving the excitable-tissue membrane potential
from its resting value
to
through a stimulation of fixed duration
satisfies:
(23)
Since from eqn. (6), :
(24)
As done in the calculus of variations let us perturb the energy-optimal time-course by the infinitesimal perturbation
, where
is an arbitrary function of time and
is an infinitesimal scalar.
(25)
From eqn. (25), the integrand in eqn. (24) becomes:
(26)
From eqns. (24) and (26), and since .
(27)
The necessary condition for to have a minimum at
for any
is:
(28)
To deal with the term of eqn. (28), it is integrated by parts :
(29)
Since the perturbation respects the boundary-value problem (BVP) with known initial and terminal conditions for
- i.e.
, then the first RHS term above vanishes. Hence, the only way that eqn. (29) will hold for any
is that we have the Euler-Lagrange-type equation:
(30)
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).
The necessary conditions for optimality require that.(32)
From eqns. (32) and (31) . Then from eqn. (33).
which is the same as eqn. (30).
From eqns. (6) and (30) we have that.and thence:
And finally, from eqn. (6).(34)
Equation (34) is a rather simple system of ordinary differential equations (ODE) that can readily be solved for a given current model to compute the energy-optimal membrane voltage profile
. The energy-efficient current waveform
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 LAP-optimal waveforms and
, 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 in eqn. (34) with
from eqn. (8):
(35)
is the membrane’s time constant and for expediency
and
= 0.
The general solution of eqn. (35) is:(36)
Given the boundary conditions and
:
(37)
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 - i.e. it has a linear rise, especially with
. Here
= 100
and
= 1
(computed using typical values from the literature for
= 1
and
= 1
).
Figure 4 presents the LAP energy-optimal stimulation profiles and
for a short and a long stimulus duration
and three membrane time constant
values.
Before we go on, it is useful to investigate the conditions for a growing exponent () waveform to outperform the
waveform.
First, has a very rapid rise. Hence, its optimal duration
will be short. Second, it is noteworthy that in [33]
= 30.4 micro-seconds! Hence, injected current rapidly leaks out. However even with the above extreme
value, at its optimal duration
the
wave does just 22% worse, which means that the
is among the best candidates for its robustly good performance.
Second, in multiple cases, the energy-optimal LAP waveform looks a lot like a ‘classical’ rectangular waveform. From eqn. (8), we may also see that, with
= 0,
= 1, the max. value of
is equal to 1 and is attained as the membrane potential reaches the threshold
. If we then replace
in eqn. (6), we see that a waveform
- that brings
from
to
at a constant rate, is the time-constant waveform
. For this example,
, which explains why
is that close to a rectangular waveform.
As a matter of fact, for very short stimulation times, the tend to be high, while
tends to be linear. Hence, the ‘classic’ rectangular (or square,
) 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 , but also of
.
Comparative properties the 
growth profiles.
The waveform may be an
waveform in disguise. I.e. some linear growth of the membrane voltage may still fit the one obtained upon ES with a
. The motivation for this is in eqn. (36), where the first term vanishes with
.
Finally, the total electric charge conveyed by the ES source may have to be considered. For example, in the of eqn. (8) the total charge consists of a capacitive charge to raise the membrane voltage by a given amount (to
), and resistive charge
. A similar situation occurs in the
model due to the opposing axial currents.
So let us solve the following auxiliary problem:
Find a linear fit to the growing exponent
, so that the ES source conveys the same resistive charge in the time interval
. I.e. we want that:
Here, for simplicity (and without any loss of generality) we have assumed and
.
For example with , we obtain
, i.e. the linear-growth equivalent has more than twice shorter duration - e.g. with
,
.
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 growth profile depend on the
ratio. As seen, for
, the optimal
is close to rectangular, while with
, the
is in effect equivalent to doing nothing for at least half of the duration, and then to a
waveform of at least doubled amplitude.
With quite similar reasoning, one can demonstrate that a 1st-order membrane voltage growth profile in the time interval
is suboptimal and equivalent to linear growth, which has about twice longer duration.
Izhikevich model.
Replacing in eqn. (34) with the
approximations from eqn. (14) or (15), see Box in Fig. 5:
(38)
The approximation of the ionic current is used for a case of very short duration (
= 10
) and the
approximation is used for a case of long duration (
= 5
). It is important to notice that - as with the
model above,
, where
(see the Box) Box: Resting-state
and asymptotic-state
ionic currents for the 0D IM; Markers are inserted at the resting and threshold membrane-voltage points, respectively
= −70,
= −55 and
= −50
.
As in the preceding model . Note that the dynamics of eqn. (38) has all FP’s of
, as well as a third FP at
, contributed by the derivative term
.
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 and
.
Figure 5 illustrates the energy-optimal LAP solution and the corresponding membrane voltage profile
. The
approximation of the ionic current is used for a case of very short duration (
= 10
) and the
approximation is used for a case of long duration (
= 5
).
It is important to notice that - as with the model above,
, where
(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:(39)where the nominal
= 1, and
.
To see how the optimal ES is affected by the level of opposing current, it is more than tempting to experiment with different values.
Hence, 3 cases are plotted in Fig. 5 - for the nominal
(cyan traces) and two additional cases: the opposing current
is either doubled (
= 2, red traces) or decreased two-fold (
= 1/2, black traces). As could be intuitively expected from the general equation (24), when
(very low ionic currents):
(40)
By the Cauchy-Schwartz inequality in the space of continuous real functions, it is straightforward to show that the voltage trajectory that minimizes eqn. (40) is such that
, where
is determined from the boundary conditions satisfied by
. Hence:
(41)
Just as in the preceding model, it is also with the shorter durations - which justifies the use of the resting approximation
.
HHM.
Here the of eqn. (34) is replaced with the resting-state -
, or asymptotic-state -
ionic current approximations (see the Box in Fig. 6).
As with the IM, bvp4c was used to numerically solve the BVP of eqn. (34). The figure follows a quite similar format to Fig. 5. 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
. The 3 solutions shown correspond to the ionic current at
(cyan trace), twice higher (thin red dash-dot), or twice lower (thick dashed black)
respectively. From eqn. (42) we can see that
= 1.6047 (half the nominal) at
, and
= 6.4188 (twice the nominal) for at
. Box: Resting-state
and asymptotic-state
ionic currents for the 0D HHM; Markers are inserted at the resting and threshold membrane-voltage points, respectively
= −77
,
= −64.55
and
= −52.35
.
Toward the gate-state variables are factored out as follows: The fast state
, while the slower variables
, and
are approximately at rest, assuming very short durations. Conversely, and assuming very long durations, toward
all gate variables are approximately at their asymptotic value, corresponding to a given membrane voltage
(see Methods).
As with the IM, we used bvp4c to numerically solve the BVP of eqn. (34) with boundary conditions and
.
Figure 6 follows a very similar format to Fig. 5.
Similarly to eqn. (39) above, 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
:
(42)where
and
= 23°C. Hence with
= 37°C, according to eqn. (42)
= 3.2094. Let this be our standard case (
= 1).
As we did with the IM, 3 gain cases are plotted in Fig. 6 for . For the two additional cases the opposing current
is either doubled (
= 2, red traces) or halved (
= 1/2, black traces).
Once again - as with the and
models above,
(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 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 functional (a function of functions), namely the energy cost (17) and the penalty (18) may conflict each other. When the penalty gain
in (18) is too low, the search will identify a lower-energy solution
, which however does not bring the membrane potential
up to the desired threshold value - i.e.
. Conversely, a too high penalty gain
will identify a very high-energy solution
, which is not only costly, but the membrane potential may also overshoot the threshold, since the ‘getting there’ is underestimated for the sake of the very last simulation steps.
As seen from Fig. 7 Panel B (which uses the approximation of the ionic current for the relatively long duration
= 2
), the linear growth profile is a reasonable estimate for the optimal membrane voltage profile
. Hence:
(43)where
is given by eqn. (41). When
is close to the LAP estimate
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
and the resulting
. With the LAP estimate, the FHOC approach resulted in a final membrane potential reasonably close to the desired threshold value - i.e.
, even if the IM was simulated with the discretized LAP waveform
(
= 10
).
panel A: discrete-time IM and FHOC panel B: continuous-time IM and FHOC, using CVODES adjoint sensitivity analysis capabilities upper plots: (dashed black) a rectangular pulse with amplitude ; (thick cyan) the LAP
; (thick black) the best FHOC
lower plots: (dashed black) linear-growth evolution of the membrane potential from
at
to
at
; (dotted gray) the desired threshold value
= −50 mV; (thick cyan) the resulting LAP
; (thick black) the resulting FHOC
.
The black traces illustrate the FHOC solution, computed for two different choices. For Panel A,
was chosen to be all zeros. When all time-step entries
were chosen to be equal to the upper bound
= 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 meta-parameter had to be kept high (
= 70) in order to respect the terminal constraint of
.
The total energy costs (all expressed as 2-norms of the obtained best ) are respectively 161, 153.2 and 423.4 (for the discrete-time version) 186.7, 159.1 and 334.2 (for the continuous-time version).
Comparing these to = 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
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 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 and
. We still need to have a very good idea of both
and
to successfully solve for
, and thence for
, in a given particular situation.
We identify also the following key and practice-oriented optimality principles resulting from the LAP perspective.
- The optimal sub-threshold membrane potential growth profile with relatively short durations
and low membrane conductivity:
First, in all simple models we used up to here, the solutionof the ODE system, defined by eqn. (34), is quite close to a linear growth from
to
. Second, with the total current
(e.g. low leak), then from eqn. (6), it follows that
will be exactly proportional to the rate of change of the membrane’s potential
. If
, then
is close to a
waveform.
- The energy-efficient waveform depends directly on the temporal shape of currents at the AP initiation site.
- The targeted
membrane voltage threshold depends on stimulation duration, with a tendency to increase with
.
- The exponential growth membrane voltage profiles
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:(44)
Without loss of generality, we used the variable 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 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
.
Hence, we may obtain the generalization of eqn. (34) as:(45)
Like the extended eqn. (44), eqn. (45) is a partial DE, depending on both temporal and spatial boundary conditions. In particular, becomes a function of
. 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 
.
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 for
ranging from the resting value (−80 mV) to some threshold
.
In this work stimulation is assumed to be intracellular and at just one spatial location (, the center RN, see Methods) along the cable structure.
To suggest a version of optimal waveforms for the MRG’02 model, we first estimate the membrane voltage threshold for each duration. One analytic way toward such estimates is readily provided by the MRG’02 model. Recall also that with simpler models
showed a tendency to increase with
.
Figure 8 presents a family of ionic current approximations at the target site (
), for a set of durations
. For each of the durations we assume that the membrane voltage trajectory
evolves according to a linear ramp from rest
to threshold
. As the latter is unknown, we produced one such ramp for each
value on the horizontal (independent-variable) axis of the figure, and then computed the corresponding ionic current
as described next.
For each of the durations it is assumed that the membrane voltage trajectory evolves according to a linear ramp from rest
to threshold
(the unknown). For each
value on the horizontal (independent-variable) axis of the figure, a
ramp was assumed and the corresponding ionic current
was computed, based on approximate gate states (see the Box). Note: for the sake of better visibility, a
gain is applied to the approx.
for the case of
= 5
. Box: For a chosen
= 5
and as
is linearly ramped up, for each gate state the plots show the ratio
, where
is given by eqn. (46) to its asymptotic value - both functions of
. Legend for gate states: opening
and closing
gates for the fast
ion-channel subtype;
persistent
channel gates;
slow
gates.
Toward gross estimates of , we first solve approximately eqn. (10) for each gate-state:
(46)where
is the gate-state value at rest and
is the average excursion from the resting membrane voltage.
Figure 8 shows the obtained approximate ionic currents as a function of just
for three very different durations -
= 0.02, 0.5 and 5
. For
= 5
, the Box in the same figure illustrates the estimated proportions-to-rest
for each of the 4 gate-state variables, at the end of stimulation.
Why does such an analysis provide upper bounds on ?
First, from the Box of Fig. 8 we can see that indeed the dynamics of the fast ion channel subtype evolves before that of the other ion channels. Particularly, we see that the estimate for inactivating
gates suggests they are completely closed for
= 5
and once
reaches around −40
.
On the other hand from the main Fig. 8, one can see that this analysis gives the intervals in which the approximate ionic currents
(i.e. remain depolarizing).
Clearly if is not reasonably within
, no miracle would yield an AP at the target location, since
becomes repolarizing outside of these bounds.
Interestingly, the analysis also predicts lowering of with longer durations. This result is exactly the opposite of what was observed with the simpler models of the HH-type, where
was repolarizing for
.
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 for the MRG’02 model. Each series was computed for the same set of 9 durations
= 20, 50, 100, 200, 400 and 500
; 1, 2 and 5
(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 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
after the end of stimulation.
With , where the minimum
was found (with 0.001 mV tolerance) using the same type of golden-section search algorithm as per the optimal
amplitude.
And the three LAP-driven series were:
The corresponding ES waveforms were computed from eqn. (44) with
.
The MRG’02 model: numerical results.
Figure 9 and table 7 illustrate the obtained as a function of
.
The computed optimal values of 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 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 optimal-waveforms’ energy and charge-transfer
values as a function of
.
As in Fig. 8 notice that the and
values are quite similar for the linear and exponential cases, for
respectively 2 and 5 ms; and also for the
-order and linear cases, for
respectively 0.2 and 0.5 ms. Toward the
values electrode impedance of 1
is assumed. Contrasted:
stands for the square (or rectangular) stimulation waveform.
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 . 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 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
= 2
, Fig. 3 (Panels A, C) and the optimal exponential growth with
= 5
, Fig. 3 (Panels B, D).
Note that we expect from an approximately globally optimal stimulation waveform to yield a specific distribution of membrane voltages
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
even when the corresponding waveform
is globally optimal. This is due for example to the small spatial constant
, which controls the spatial diffusion with time.
However, if the spatial profile is about the same for different durations and the corresponding different waveforms
(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
= 5
have duration
= 2.3
. 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 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 re-polarizing node currents for each different duration. The waveforms’ offsets (due to ) are inversely proportional to duration. This readily compares qualitatively with the results in [8]. Especially for very short durations (e.g.
), the optimal waveform
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 (
).
Since here , where
is given by eqn. (41), from eqn. (6)
. The figure is optimized to present clearly both
and
(*1) The dashed trace at the bottom plots
as a function of
(*2) Toward equally good plot visibility, for all durations
, the waveforms
are rubber-banded to take the same graph width as the 1 ms-waveform. This is illustrated by the scale bars for the shortest duration
= 20
. (*3) The vertical scale is the same for all plots, except for the logarithmic offset, as defined by pt. (*1) above.
Numerous essential differences in the approach preclude further objective comparisons. Interestingly however, for the longer durations ( 0.5
) the results in [8] show very little (if any) variation with
(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 from its resting (
) to its threshold value (
) through a stimulation of fixed duration. Through the LAP perspective, we obtained eqn. (34) - a general (model-independent) description of the energy-optimal time-course of the excitable-tissue’s membrane potential
.
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 . 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
), 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
and hence there is no danger of
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 optimization 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 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 values of the optimal rectangular waveform, for
= 100 and 200
- see Table 9. As can be seen from Fig. 9, no benefit in terms of lower
can be associated to the steep rise of the rectangular waveform, since
is expected to be higher, esp. for dramatically shorter durations. This was further confirmed by numerical testing with dual linear (high/low rate)
rise schedules (data not shown), which all had inferior performance to the baseline simple linear-growth protocol. However, the rectangular waveform also leads to steep capacitive decay of
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 the ones associated to the linear-ramp voltage profile (see Fig. 10). On Fig. 13 one can see that the steep rise of the
waveform yields an early super-linear ramping of the membrane voltage. However, the rectangular waveform requires a lot more charge
to be transferred.
One can see clearly that the steep rise of the waveform yields an early superlinear ramping of the membrane voltage. However, the rectangular waveform requires a lot more charge
to be transferred (see Fig. 10).
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.
Author Contributions
Conceived and designed the experiments: NIK SMD AV. Performed the experiments: NIK. Analyzed the data: NIK SMD AV FR. Contributed reagents/materials/analysis tools: NIK SMD AV. Wrote the paper: NIK SMD AV FR MS.
References
- 1. Liberson W, Holmquest M, Scot D, Dow M (1961) Functional electrotherapy: stimulation of the peroneal nerve synchronized with the swing phase of gait of hemiplegic patients. Arch Phys Med Rehabil 42: 101–105.
- 2. Dhillon G, Horch K (2005) Direct neural sensory feedback and control of a prosthetic arm. IEEE Trans Neural Syst Rehab Eng 13(4): 468–72.
- 3. Fetz EE (2007) Volitional control of neural activity: implications for brain-computer interfaces. Journal of Physiology-London 579: 571–579.
- 4. Doherty J, Lebedev M, Hanson T, Fitzsimmons N, Nicolelis M (2009) A brain-machine interface instructed by direct intracortical microstimulation. Front Integr Neurosci 3: 20.
- 5. Lucas T, Fetz E (2013) Myo-cortical crossed feedback reorganizes primate motor cortex output. J Neurosci 33(12): 5261–74.
- 6. Sahin M, Tie Y (2007) Non-rectangular waveforms for neural stimulation with practical electrodes. J Neural Eng 4(3): 227–33.
- 7. Foutz T, McIntyre CC (2010) Evaluation of novel stimulus waveforms for deep brain stimulation. Journal of Neural Engineering 7: 066008.
- 8. Wongsarnpigoon A, Grill WM (2010) Energy-efficient waveform shapes for neural stimulation revealed with a genetic algorithm. Journal of Neural Engineering 7(4): 046009.
- 9. Wongsarnpigoon A, Woock JP, Grill WM (2010) Efficiency analysis of waveform shape for electrical excitation of nerve fibers. IEEE Trans Neural Syst Rehab Eng 18: 319–328.
- 10. Weiss G (1901) Sur la possibilite de rendre comparable entre eux les appareils servant a l’excitation electrique. ArchItalBiol 35: 413–446.
- 11. Lapicque L (1907) Recherches quantitatives sur l’excitation electrique des nerfs traites comme une polarization. J Physiol (Paris) 9: 622.
- 12.
Lapicque L (1926) L’excitabilite en Fonction du Temps. La Chronaxie, sa Signification et sa Mesure. Presses Universitaries de France, Paris.
- 13. Lapicque L (1931) Has the muscular substance a longer chronaxie than the nervous substance? The Journal of Physiology 73: 189–214.
- 14.
Lykknes A, Opitz DL, Tiggelen BV, editors (2012) For Better or For Worse? Collaborative Couples in the Sciences, volume 44 of Science Networks. Historical Studies. Basel: Springer-Birkhaeuser, 319 pp.
- 15.
Rattay F (1990) Electrical Nerve Stimulation: Theory, Experiments and Applications. Wien, NY: Springer-Verlag.
- 16. Geddes LA (2004) Accuracy limitations of chronaxie values. IEEE Trans Biomed Eng 51: 176–181.
- 17.
Feynman R, Leighton R, Sands M (1964) The principle of least action. IN: The Feynman Lectures on Physics Addison-Wesley, vol.II, ch.19.
- 18.
Rall W (1977) Core conductor theory and cable properties of neurons. In: Handbook of Physiology - The Nervous System (APS) vol 1, chap 3, 39–97.
- 19. Katz B (1937) Experimental evidence for a non-conducted response of nerve to subthreshold stimulation. ProcRoySocLondon serB 124: 244–276.
- 20. Rattay F, Wenger C (2010) Which elements of the mammalian central nervous system are excited by low current stimulation with microelectrodes? Neuroscience 170: 399–407.
- 21. Hu W, Tian C, Li T, Yang M, Hou H, et al. (2009) Distinct contributions of Nav1:6 and Nav1:2 in action potential initiation and backpropagation. Nat Neurosci 12: 996–1002.
- 22.
Izhikevich E (2007) Dynamical Systems in Neuroscience. Cambridge, MA: MIT Press.
- 23. Fitzhugh R (1961) Impulses and physiological states in theoretical models of nerve membrane. Biophys J 1: 445–66.
- 24.
Izhikevich E (2004) Which model to use for cortical spiking neurons? IEEE Transactions on Neural Networks 15(5).
- 25. McIntyre CC, Richardson AG, Grill WM (2002) Modeling the excitability of mammalian nerve fibers: Inuence of afterpotentials on the recovery cycle. Journal of Neurophysiology 87: 995–1006.
- 26. Ranck JJB (1975) Which elements are excited in electrical stimulation of mammalian central nervous system: A review. Brain Research 98: 417–440.
- 27. Rattay F (1998) Analysis of the electrical excitation of cns neurons. IEEE Trans Biomed Eng 45: 766–772.
- 28. Rattay F (1999) The basic mechanism for the electrical stimulation of the nervous system. Neuroscience 89: 335–46.
- 29.
Danner S (2010) Master thesis: Computer simulation of electrically stimulated nerve fibers in the human spinal cord. TU Vienna, Inst for Scientific Computing & Analysis, Advisor: F Rattay.
- 30. Danner SM, Hofstoetter US, Ladenbauer J, Rattay F, Minassian K (2011) Can the human lumbar posterior columns be stimulated by transcutaneous spinal cord stimulation? A modeling study. Artificial Organs 35: 257–262.
- 31.
Grill W (2010). Personal comm. w.r.t. possible printing errors in the MRG’02 model’s original 2002 publication.
- 32.
Sage AP, Melsa JL (1971) System Identification, volume 80 of Mathematics in Science and Engi-neering. Academic Press, 221 pp.
- 33. Jezernik S, Morari M (2005) Energy-optimal electrical excitation of nerve fibers. IEEE Trans Biomed Eng 52(4): 740–3.
- 34. Shaw RM, Rudy Y (1997) Ionic mechanisms of propagation in cardiac tissue. roles of the sodium and l-type calcium currents during reduced excitability and decreased gap junction coupling. Circ Res 81: 727–41.
- 35. Rattay F, Paredes LP, Leao RN (2012) Strength-duration relationship for intra- versus extracellular stimulation with microelectrodes. Neuroscience 214: 1–13.
- 36. Hofmann L, Ebert M, Tass PA, Hauptmann C (2011) Modified pulse shapes for effective neural stimulation. Frontiers Neuroeng 4: 9.