Fundamental Limits of Forced Asynchronous Spiking with Integrate and Fire Dynamics

Received: 30 January 2017 / Accepted: 25 September 2017 / © The Author(s) 2017. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.


Introduction
The manipulation of networks of neurons in the brain through the use of extrinsic controls-neurocontrol-is a key problem in experimental neuroscience [1]. Such capability has the potential to enable new and important study of questions in neural coding or how the firing activity of brain cells determines their ability to carry and process information [2]. Moreover, improving the use of neurostimulation may aid the refinement of how such technology is used in clinical settings [3,4].
The use of stimulation in the study of neural coding is itself an established paradigm in neuroscience. The general idea is straightforward: by inducing neural activity and observing the consequent behavior of the organism, we can infer the functional role of the region in question. For example, cortical microstimulation of certain brain regions has been shown to induce behavioral changes in the context of perceptual tasks such as visual decision-making [5,6]. Recently, several key advances in neurostimulation technology, such as the advent of optogenetics [7], have made neurocontrol possible at unprecedented spatial scales. Thus, experimentalists are able to assess the functional role not simply of different neural populations, but potentially of specific neurons and the timing of their spikes. That is, it may be possible to test the long-standing neural coding hypothesis that spike timing is crucial to information processing [8].
Currently, however, these hardware instantiations are typically used in perturbative paradigms wherein "pulses" of input are used to alter neural firing in a bulk manner (see Fig. 1) that does not control the precise timing of individual neuronal spikes. Formal control analysis or design in this context, though desired, is not well studied [9]. Thus, there is a need for formal mathematical analysis regarding the fundamental limits of such stimulation, particularly as it pertains to the feasibility of inducing precisely timed spiking activity in neural populations (Fig. 1).

Prior Work in Neuronal Control
The control of neural activity has received substantial attention in the context of oscillations and synchronization, spurred in large part by interest in clinical brain stimulation for motor disorders [10,11]. The objective in this class of neurocontrol problem is generally the forced splaying of neural phases (i.e., desynchronization), The use of such stimulation has historically been limited to perturbative paradigms, wherein pulse-type inputs are used to create bulk population responses without fine temporal structure. (B) Increasingly, experimentalists seek to induce more precise spiking patterns in specific subsets of the population, which may necessitate the design of nuanced stimulation waveforms.
In contrast, we consider herein the mathematical problem of asynchronous neurocontrol (i.e., control neural spiking without overt rhythmicity), iIn other words, forcing a neuron to spike but not necessarily periodically. The other key distinction of our work is that we consider a neuronal-level objective (i.e., spiking and spike timing) versus a population-level objective (i.e., synchronization or desynchronization). We have previously provided early formulations of this problem and highlighted key analytical challenges in the development of controllability analysis for spiking models [23,24]. Other works regarding formal control design include optimal control design for a single neuron [25] and using statistical modeling frameworks [26,27].

Neurocontrol with Common Input
A key challenge associated with neurocontrol is underactuation, wherein a small number of inputs (in many current implementations, a single input) impinges on an orders-of-magnitude greater number of neurons [23], as schematized in Fig. 1. In other words, individual neurons are not addressed via independent inputs, but rather a common one. This challenge is ubiquitous across stimulation modalities and is, perhaps, the major constraint that has restricted the use of neurostimulation to the aforementioned perturbative paradigms. In the context of the discussed oscillatory objectives, some progress has been made on solving control problems such as entrainment and synchronization in the presence of underactuation [28][29][30][31]. However, this issue is unresolved in the case of asynchronous timed spike control objectives, such as those we consider herein. Current and foreseeable neurostimulation technologies are likely to face the challenge of underactuation, especially for in vivo instantiations.

Specific Contributions
In this paper, we address the problem of time-optimal control of spiking in pairs of Leaky Integrate-and-Fire (LIF) neurons, where the desired spiking is selective, that is, certain neurons spike while others remain silent. We specifically focus on the case where two neurons receive a common input, which, as mentioned before, is a key constraint in the practical design of neurocontrol methods. Our major contributions are in the characterization of fundamental limitations for neuron-level control as revealed through a formal mathematical analysis. This treatment leads to the postulation of practical neurocontrol design strategies. Specifically, we provide: 1. The formal synthesis of time-optimal selective spiking solutions in pairs of LIF neurons. The synthesis involves application of the Pontryagin maximum principle, but with several nontrivial caveats due to the selectivity specification, which leads to state constraints. We prove that the optimal solution in this case involves use of the so-called boundary control, associated with the state constraints. Sufficient conditions for optimality are verified.
2. The formal synthesis for time-optimal control of longer sequences of spikes. Here, the solution is derived via dynamic programming, but again with several nontrivial developments due to nondifferentiability of the value function. In particular, we prove the nonexistence of an optimal solution for specific classes of spike sequences. 3. The development of design methods for timed patterns of spikes. In this case, there is no unique optimal solution. Nevertheless, we derive a greedy algorithm that can provide near-perfect construction of patterns under specified conditions. Finally, we evaluate the performance of our control design when the system is subjected to noise and disturbances.
Our presentation and discussion on fundamental optimal control analysis and design work toward the overall goal of understanding the limits of neurocontrol. We illustrate several interesting control phenomena that arise due to the peculiarity of spiking dynamics. Specifically, the problem considered, although ostensibly simple, leads to several interesting features in the optimal control synthesis due to state constraints.

Definitions: Spike Sequence and Pattern Control
We begin by formally defining the notions of spike sequences and patterns, which will facilitate our approach to spike timing control.
Definition 1 (Spike Sequence) In a population of N neurons, an M-spike sequence is a vector where σ k ∈ {1, 2, . . . , N} indicates the neuron that produces the kth spike in the sequence.
Definition 2 (Spike Pattern) In a population of N neurons, an M-spike pattern is a sequence with timing, that is, where σ k ∈ {1, 2, . . . , N} indicates the neuron which produces the kth spike at time The goal of this paper is to provide a set of fundamental characterizations regarding the time-optimal control of spike sequences and patterns.

Model Formulation
We proceed with the model formulation, starting with the base model and then adding synaptic coupling between neurons. Fig. 2 The leaky integrate-and-fire circuit. The membrane potential rises under the stimulus u(t) until it hits the threshold V T , at this point v(t) is artificially reset to V rest and a spike is said to be generated. We also show a cartoon of the possible voltage trace of the neuron under a rectangular pulse input.

Base Model
The integrate-and-fire neuron is a well-established model in computational neuroscience [32,33]. The circuit of this model is shown in Fig. 2, where a capacitor C and resistance R (modeling the capacitive and resistive properties of the cell membrane) are in parallel, with u(t) being the external stimulus. Denoting the membrane potential as v(t), the charge deposited on the capacitor is q = Cv, and therefore the current is given by I C = C dv dt , leading to the linear dynamics where V rest is the resting potential, and κ m = RC is the membrane time constant.
Here, I syn denotes synaptic input entering from other neurons. We also introduce a parameter β that encapsulates the effectiveness of the external input u(t) for each neuron. Spike generation. In this model, a spike is said to be generated at time t s if the membrane potential reaches a predetermined threshold voltage V T . Upon emitting a spike, the membrane potential is reset to V rest . Thus, spike generation is governed by the discontinuous resetting rule Model normalization. In what follows, we assume that V rest = 0. This normalizing assumption is not restrictive, since it can be readily achieved by a simple translation in the coordinate system, that

Synaptic Coupling
We build an approximate model of synaptic coupling based on the standard formulations in [33]. Key to this formulation is the notion of impulsive coupling, wherein the major effect of I syn occurs during a brief time window following an afferent spike (i.e., a spike from another neuron). Following a reduction of continuous synaptic models (see Appendix A.1), we formulate I syn as where T denotes the set of all afferent spike times, and ρ syn (t) is a synaptic constant that depends on the specific parameters of the neuron. If all neurons remain below the threshold, then I syn ≡ 0. Thus, the effect of a synaptic event on the postsynaptic neuron can be understood as an instantaneous rise in voltage that occurs only when a neighboring, connected neuron fires a spike. Knowing this rise can allow us to insulate neurons from each other in the spike control problem, formulated in the next section.

Problem Formulation: Minimum Time Selective Spiking
In this paper, we study three base problems pertaining to the design of u(t) to create structured spiking patterns in populations of two LIF neurons of the form (3). We first consider the problem of time-optimal sequence control, that is, inducing target sequences with minimal temporal spacing between the beginning and end of the sequence. It turns out that this problem amounts to an analysis of selective spiking. We formulate a canonical version of this problem in two dimensions.
Problem 1 (P1: Pairwise time-optimal selective spiking with synaptic guard) Consider two coupled LIF neurons of the form (3): where and I syn i are impulsive synaptic inputs of the form (5) for i = 1, 2. Find the control input u(t) such that with arbitrary initial condition v(0) ∈ G, where and u(t) solves the time-optimization minimize J(u) = τ 0 dt (9) over all measurable functions u that take values in the control set, where U is this set of admissible inputs.
Functional decoupling of the network via guard V G . The parameter V G in (7), referred to as a synaptic guard, is key to selectivity. It ensures that Neuron 2 remains below threshold and, further, is insulated from the synaptic effect due to the induced spike in Neuron 1, that is, where ρ syn (t) is the synaptic contribution to the post-synaptic neuron (here, Neuron 2) and is derived in Appendix A.1. The guard, in essence, keeps the nonselected neuron sufficiently away from its own threshold so as not to produce an undesired, collateral spike. It is important to note that in solving (P1), it is sufficient to consider the dynamics in (6) asv since both neurons are below threshold for the duration of the synthesis. Despite this simplification in the dynamics, the selectivity/guard criterion (7) poses a key challenge, that is, it is not sufficient to simply fire Neuron 1 in minimum time, since doing so may in general cause Neuron 2 to fire an undesired spike. Mathematically, (7) functions as a state constraint that, as we will see, leads to several complications in the optimal synthesis. If the problem has a solution for either choice of neuron labeling, then the population is said to be pairwise feasible, that is, either neuron can be made to spike selectively.
Problem 2 (P2: Pairwise time-optimal selective sequencing) For the two-neuron network in (11), find the control input that achieves any M-spike target spike sequence Σ S time optimally, that is, The key complication here is the nondifferentiability of the value function within the dynamic programming, as well as the spike discontinuity (4).
Problem 3 (P3: Pairwise time-optimal selective patterning) Considering the same model in (11), find the control that induces the spiking in the two neurons according to the times specified in the target pattern Σ P , constrained by the underlying sequence. Mathematically, with the same constraints as described in (13) and t 0 = τ 0 = 0. Note that t k are the desired spike times, and τ k are the actual spike times.

Minimum Time Selective Spiking
We consider the minimum-time selective spiking problem P1. We assume, without loss of generality, that the neurons are labeled so that the objective is to fire Neuron 1. It turns out that the solution to this problem depends on the ratio (see Appendix A.2) which we treat in two separate cases corresponding to γ 1 ≶ V T V G . As we will show in the following sections, for γ 1 > V T V G , that is, Case 1, selective spiking can always be accomplished. However, if γ 1 ≤ V T V G , that is, Case 2, a solution may not exist, and pairwise feasibility is not guaranteed. (11), where

Proposition 1 Consider the two-neuron network
Assume that the set of admissible controls U forms a box constraint of the form U = [0, U], and we take as given the initial conditions v i (0) < V G , i = 1, 2. The time optimal feedback control u * ∈ U for the selective spiking problem P1 for Neuron 1 is given by where u arc = a 2 b 2 V G is the unique control that keeps v 2 (t) = V G invariant. Moreover, such a control always exists. Thus, optimal controls are either given by a constant control at maximum value, u * (t) ≡ U , if the state space constraint does not become active, or if the corresponding trajectory meets the state space constraint, then optimal controls are a concatenation of a segment for the maximum control until the state constraint is reached followed by a constant boundary control u * (t) = u arc until the terminal value v 1 = V T is reached.
Proof Necessary conditions for optimality for problem P1 are given by the Pontryagin maximum principle. In the presence of state space constraints, these take a rather complicated form (the multipliers associated with the state space constraint are measures). The problem considered here, however, is simpler, and instead of analyzing those conditions, we shall define a synthesis of extremal controlled trajectories through a direct construction and then verify the optimality of the synthesis. In particular, there is no need to consider possible degeneracies that in principle are allowed by necessary conditions for optimality (e.g., abnormal extremals, etc.).
Synthesis Construction. We want to solve the optimal control problem P1 on the set G in (8). We first treat the problem in the absence of the state constraint and define the Hamiltonian function as According to the maximum principle, as long as no state space constraints are active, the multiplier λ is a solution to the adjoint equatioṅ and the optimal control minimizes the Hamiltonian over the control set [0, U]. The solutions of (19) are of the form for some constants c 1 and c 2 , and thus with as the switching function. The terminal constraint is defined by ψ(τ, v) = v 1 (τ )−V T , and the transversality condition [34, Sect. 2.2] of the maximum principle implies that λ(τ ) = [ν 0] where ν is some multiplier. This gives us c 2 = 0, and thus the switching function has a constant sign in the absence of the guard constraint. Hence the optimal control is simply a BANG, that is, the maximal input. With the state constraint (the guard), there can be switching in the optimal control, and we need to consider two subcases: trajectories that do or do not hit the boundary v 2 = V G . For A with real eigenvalues, the optimal controls of linear single input control systems are BANG-BANG with at most n − 1 switchings (where n is the dimension of the system; here n = 2) [34], and we must have u > 0 at the spike time (otherwise, v would be decaying). We thus consider controls only of the form These define a smooth flow of extremal controlled trajectories as long as the state space constraint is not violated. If the extremals hit the state constraint boundary, then the control must switch to the boundary control u arc that keeps the system from exceeding the constraint: However, we need to verify whether this boundary control u arc will eventually bring Neuron 1 to threshold. For v 1 = V T and u = u arc , we havė where the inequality holds by our assumption on γ 1 . Now, if (25) holds, then in facṫ v 1 > 0 for all v 1 ∈ [0, V T ] under the boundary control, and v 1 will eventually reach threshold. Thus for appropriate initial conditions, applying the maximal input u(t) = U produces a spike in Neuron 1 without hitting the Neuron 2 guard. For the remaining initial conditions, we construct a control that applies maximal input until the guard is reached and then drops to u arc until v 1 hits threshold. Note that we do not need to employ the zero control in (23), so we may taket = 0 (the possibility of additional switching will arise in the next section under the alternative case for γ 1 ). Thus the control (17) will produce a spike in Neuron 1 without inducing a spike in Neuron 2 across all initial conditions. This concludes the synthesis construction.
Proof of Optimality. The optimality of this control follows from regular synthesistype sufficient conditions for optimality, and we briefly outline the reasoning. The value or cost-to-go function of this synthesis is continuous but not differentiable on the curve that separates initial states for which the trajectory includes a boundary segment from those that do not. The curve Γ that separates these two regions is defined by the set of initial conditions that hit the final condition To find this curve, we first explicitly compute the time for v 1 to hit threshold, where for convenience we define We then eliminate τ by solving explicitly for v 2 to find the separatrix as We define the region Γ − as bounded between Γ and v 1 = V T inclusive, and the region Γ + = G \ Γ − . Thus, Γ + includes all initial conditions whose trajectories include a boundary arc, whereas initial conditions in Γ − can be driven to threshold directly at maximum input.
The value function corresponding to this synthesis is For trajectories without a boundary arc, the value is just the spike time under maximal input, calculated as in (26), The calculation of the value V + (v) involves two steps: the time t g for Neuron 2 to reach the guard voltage, plus the time t th for Neuron 1 to attain the threshold V T under the boundary arc control. By direct calculation, where is the Neuron 1 voltage at the time t g , that is, when the trajectory hits the Neuron 2 guard. It is clear from the construction that V is continuously differentiable in the interior of G away from the curve Γ . We now show that on Γ , V remains continuous, but is no longer differentiable. Substituting v 2 from (29) into (33) yields Hence (32) reduces to Substituting v 2 once again into (35), it follows that However, so that V is not continuously differentiable. All controlled trajectories in the synthesis are extremals, and away from Γ , the value function V satisfies the Hamilton-Jacobi-Bellman equation for the unconstrained optimal control problem where L is the Lagrangian of the problem (for time optimal control problems, as in our case, L = 1). This conclusion follows from the method of characteristics (e.g., see [34]) but can also directly be verified using the explicit formulas derived above. That V is not differentiable on Γ does not invalidate the proof of optimality, although the standard optimality argument based on dynamic programming (e.g., [34], Theorem 5.2.1) does not apply. Here, we need to invoke regular synthesis constructions (see Appendix A.3) as they are described in [34,Sect. 6.3]. Since trajectories do not return from the state space constraint into the interior of the state space, these arguments could, for example, be undertaken by redefining the state space constraint as a second terminal manifold, along with a penalty term that gives the time along the boundary control until v 1 = V T . Alternatively, the constructions in [35], where a regular synthesis argument has been generalized to problems with order 1 state space constraints, could be modified to apply to cases where the state space constraint is active at the terminal time. Either way, straightforward modifications of regular synthesis type arguments give the optimality of the above field of extremals.

Example 1
We demonstrate minimum spike time control in an example of (11) with the following parameters: Note that these are idealized parameters used for illustrative purposes only, although with biologically plausible units. Here, the condition γ 1 > V T V G is satisfied, and we can apply the above proposition to induce a spike in Neuron 1 in minimal time.

Selective Spiking, Case 2: γ 1 ≤ V T V G
We now consider the case of eliciting a spike in Neuron 1 when γ 1 ≤ V T V G . We showed in the previous section that for Case 1, a control solution always exists. It will turn out that not all parameters allow a solution in Case 2, so this case reveals the conditions for pairwise feasibility of sequences while providing the minimum time spiking solution when it exists. We might expect the solution in Case 2 to be qualitatively similar to Case 1, but in fact there are no longer increasing trajectories that ride along the guard boundary: under the boundary control (u arc = a 2 V G b 2 ), we findv 1 < 0 at v 1 = V T , that is, along the guard, v 1 (t) does not rise beyond a certain limit and fails to reach the threshold V T . Instead, we have the following: Assume that the set of admissible controls is a box constraint U = [0, U]. The time optimal control u * ∈ U for the selective spiking problem P1 for Neuron 1, if such a solution exists, is with Γ ± defined as before.
Proof We follow a similar analysis to the previous case, but identify the differences in the optimal control structure from the solution in Sect. 3.1. Again, our approach is to define a synthesis of extremal controlled trajectories, prove their optimality, and finally give conditions for the existence of a solution for all v ∈ G. Synthesis Construction. The Hamiltonian and multiplier are similar to (18) and (20). The minimum condition similarly results in (21) with the conclusion that the optimal control is simply BANG at u * (t) = U for trajectories that do not hit the guard under this control. Similarly to (28), there again exists a curve Γ that separates such initial conditions from those requiring switching, given by (29). Note that there is no boundary segment in this case as u arc cannot drive the voltage of Neuron 1 up to threshold along the state constraint boundary (see Appendix A.2), and thus we are led to consider controls only of the form in the interior of G, andt = 0 is allowed. This concludes the synthesis construction. Proof of Optimality. The value function for the region Γ − equals the time taken by Neuron 1 to reach the threshold V T under the constant control U and takes the same form as (31). For v ∈ Γ + , the value function is calculated assuming that the control is turned off for an interval [0,t], during which the system decays from the initial wherê Here we cannot get an explicit expression for V + in terms of the initial condition [v 1 v 2 ] T because of the transcendental form of (43). Note that, for this synthesis, the state space constraint does not become active. It is clear from the construction that the corresponding values satisfy the Hamilton-Jacobi-Bellman equation away from Γ . However, this problem is nonstandard in that the value function may no longer be continuous on Γ , with the only exception at v 1 = 0, that is, In general, there may exist a unique point on the curve Γ (in our problem with u = 0) where the vector fieldv = Av is tangent to Γ while pointing in the opposite direction. As a result,v = Av points into the region Γ + and into the region Γ − , above and below this point, respectively. This generates a loss of small-time local controllability that causes the value function to become discontinuous along Γ above this point. For, if the initial condition lies to the right of Γ + above this point, then optimal trajectories must decay below the point in order to reach the terminal manifold. We see this in Fig. 3(b), where the OFF segment in the extremal cannot simply converge to the separatrix Γ , no matter how close it is to Γ . This issue of controllability makes the value function discontinuous. The value is still lower semicontinuous on the full state space. In fact, the value of this synthesis satisfies Sussmann's weak continuity requirement [34,Definition 6.3.3]. Although the discontinuity of the value impedes on the application of most HJB-type sufficient conditions for optimality, this is not the case for regular synthesis-type constructions (see Appendix A.3), and the optimality of the synthesis follows from Theorem 6.3.3 in [34].
Existence of Solution. However, the control approach in (41) will fail if trajectories starting in Γ + do not in fact hit the separatrix at some time during the initial offcontrol. A necessary and sufficient condition for trajectories to hit the separatrix is that Γ intersects the positive v 2 axis. When this condition holds and v(0) lies above Γ , then there must be a timet where the trajectory hits Γ under u = 0. Conversely, suppose Γ does not intersect the positive v 2 axis. The slope of Γ , considering v 2 as a function of v 1 , must be less than the slope of the decaying trajectory for there to be an intersection (ignoring the degenerate parameter choice for which tangency is possible). Taking the ratiosv 2 /v 1 for u = 0 and u = U (recalling that Γ is itself a solution with maximal input) and rearranging the result show that the slope condition can be met only if v 2 > γ 1 v 1 . However, by our assumption γ 1 ≤ V T /V G , no point on Γ meets this inequality (the curve lies entirely below the line from the origin to [V T V G ] T ). In fact, sincev i , i = 1, 2, is monotonic in u, it follows that there is no admissible control that can push a solution across Γ , so that the latter serves as a barrier to Neuron 1's threshold for all initial conditions in Γ + (at least, without first crossing the Neuron 2 guard). So in this case, selective spiking of Neuron 1 is not possible.
Thus, the condition for the existence of a time-optimal solution for selective spiking of Neuron 1 is that the v 2 intercept of Γ is positive, which occurs when Example 2 We use the same parameter values as in (39) but swap the roles of Neuron 1 and Neuron 2, that is, Now, γ 1 ≤ V T /V G . Moreover, condition (45) holds, so that the switching separatrix intersects the positive v 2 axis. Thus a time-optimal solution for selectively spiking Neuron 1 always exists. Figure 3(b) shows example trajectories.

Geometric Interpretation of Cases and Pairwise Feasibility
Thus far in our discussion we assume, without loss of generality, that a selective spike is desired in Neuron 1. Now for pairwise feasibility, that is, to analyze when timeoptimal selective spiking of either neuron is possible (from any initial condition), both neurons must be associated with either Case 1 or Case 2. To do this, we introduce We associate Neuron 1 with γ 1 and Neuron 2 with γ 2 to determine the case (Sects. 3.1 and 3.2) to which these neurons belong. We say Neuron 1 is Case 1 or 2 when γ 1 > V T V G or γ 1 ≤ V T V G , respectively, and similarly for Neuron 2 with the same inequality relation on γ 2 . Since we have V T > V G , this allows for three possible scenarios, Neuron 1 is Case 1, and with γ 2 being the reciprocal of γ 1 , we have Neuron 2 is Case 2. 2. γ 1 < V T V G , γ 2 > V T V G : Neuron 1 is Case 2 and Neuron 2 is Case 1, and the structure of the solution is identical to the previous scenario.
Neurons are Case 2, and this happens when V G V T ≤ γ 1,2 ≤ V T V G . As we will show in the following sections, for one of the neurons belonging to Case 1, pairwise selective spiking can be accomplished. However, if γ 1,2 ≤ V T V G , that is, both neurons are Case 2, a solution may not exist, and pairwise feasibility is not guaranteed.
To provide an additional geometric interpretation (see Appendix A.2) of these conditions, we introduce the quasi-static equilibrium line which defines the set of points for whichv = 0 (for each u ∈ U ). In a pair of neurons, the following two possible parameterization scenarios can be encountered.

Neuron 1 and 2 Correspond to Different Cases
Here we discuss the pairwise feasibility for when Neuron 1 is Case 1 and Neuron 2 is Case 2. It is important to note that the result extends to the reverse scenario, that is, Neuron 1 is Case 2 and Neuron 2 is Case 1.
Here, the line of quasi-static equilibrium in (48) intersects the line v 1 = V T before it intersects v 2 = V G . Thus, Neuron 1 can always increase along the Neuron 2 guard boundary. Conversely, Neuron 2 cannot increase along the Neuron 1 guard beyond the point of intersection between v(∞) and v 1 = V G . As we showed before, in this case, selective spiking of Neuron 1 is always possible. Thus, pairwise feasibility reduces to condition (45) modulo a swapping of labels. Specifically, we have the following: Proof The proof follows immediately from Proposition 2 and (45), with a swapping of labels.
Thus, it follows that if (49) does not hold, a time-optimal solution for Neuron 2 does not exist (for all initial conditions), and thus the neurons are not pairwise feasible.

Neuron 1 is Case 2; Neuron 2 is Case 2
If both neurons are Case 2, then pairwise feasibility would necessitate (49) holding to within a swapping of labels (i.e., so that either neuron can be selectively spiked). Clearly, this is impossible (see Appendix A.2) except for the limiting case where V G = V T , that is, the neurons are not guarded. In such a scenario, the optimal solution may produce simultaneous spiking of both neurons depending on the initial condition.

Minimum Time Sequence Control
We now use the above results to analyze longer pairwise spiking sequences Σ S to solve the problem P2. Based on the results of the previous section for pairwise feasibility, that is, to allow all possible spike sequences for two neurons, we make the following assumption hereon.

Assumption 1
The pair of neurons are parameterized so that Neuron 1 satisfies Case 1, Neuron 2 satisfies Case 2, and Lemma 1 holds.
This assumption ensures that the selective spiking solutions for the two neurons are given by Proposition 1 and 2, respectively. We now analyze all the possible length 2 sequences, that is, [1,1], [1,2], [2,1], and [2,2], and recognize how we can use the basic characterizations developed in Sect. 3.1 and 3.2 to synthesize a time-optimal strategy for these sequences. We employ a dynamic programming approach where, using the time-optimal solution for the second spike in neuron i, we define a terminal cost and then solve the resulting optimal control problem for the first spike in neuron j , i, j ∈ {1, 2}. Whereas the optimal synthesis for some of these sequences can be generalized from the solution of P1, we shall see that for the target sequence [2,1], no time-optimal control solution may exist.

Synthesis of All-2 Spike Sequences
Without loss of generality, consider the spike sequence Σ S = [1, 1] that we want to achieve in minimum time. We will use the concept of dynamic programming to solve the following problem: Optimal trajectories for sequence [1,2]. (c) Optimal trajectories for sequence [2,2]. In this case, all initial conditions collapse onto a single manifold associated with the second spike.
We will start from the last spike, Neuron 1, for this example and solve the minimum time problem P1 for all the initial condition for Neuron 2, namely v 2 ∈ [0, V G ], v 1 = 0, and use the solution of P1 as the terminal cost ϕ(v(τ 1 )) for the previous spike, Neuron 1 again, in our case. So we will solve the following optimal control problem: Now we will seek synthesis for all possible two spike sequences using (51).

Spike Sequence [1, 1]
The optimal synthesis for the sequence Σ S = [1, 1] is given in Fig. 4(a). We highlight the solution of P1 for Neuron 1 on the top left, the terminal cost ϕ(v 2 (τ 1 )) in the middle, and in the bottom, we show the solution of (51). On the right, we construct the complete synthesis for the whole sequence.
Given an arbitrary initial condition [v 1 v 2 ] T , the time-optimal solution of the first part without any terminal cost (i.e., ϕ(v 2 (τ 1 )) ≡ 0, given by Proposition 1) has the property that, among all admissible controls, it leads to the smallest possible value for the terminal state v 2 (τ 1 ). Since the function ϕ(v 2 (τ 1 )) is strictly increasing, this is then also the optimal solution for the combined problem and thus allows us to simply concatenate two solutions of P1 for Neuron 1. Overall, the optimal control is simply given by the BANG control U until v 2 reaches the guard, after which the boundary control is used exactly as in the single spike problem. [1,2] However, such monotonicity arguments do not work in the other cases. Figure 4(b) shows the synthesis of optimal controlled trajectories for the sequence Σ S = [1,2]. The terminal cost ϕ(v 2 (τ 1 )) is calculated as the value function from the solution of P1 for Neuron 2 and is a strictly decreasing function of v 2 (since the higher the voltage v 2 , the lower the time to induce a spike in Neuron 2). Thus, in principle, it might be possible for the solution of the first part to deviate from the solution of P1 for Neuron 1 if the loss in doing so would be made up by the gain in the penalty function ϕ(v 2 (τ 1 )) at the terminal point. Consider the switching function

Spike Sequence
If there is a switching at t =t, then we have Also, for a switching structure OFF-BANG, we must havė Now we use (53) for computing the derivative of the switching functioṅ since the terminal cost is a decreasing function of v 2 . Also, we have previously derived that the adjoint variables are solutions of linear homogeneous differential equations that do not change sign in t ∈ [0, τ 1 ]. So we have λ 2 (t) < 0 as well. Using these and assuming that a 2 < a 1 , from (55) we geṫ This violates the necessary condition in (54) for an OFF-BANG switching. Note that for the case a 1 < a 2 , OFF-BANG switching cannot be ruled out using this argument, and the synthesis has to be constructed by direct computation. In our example with the parameters from (39), it turns out that the optimal solution is simply BANG/BANG-BOUNDARY (17), that is, the terminal cost ϕ(v 2 ) has no effect on the solution of (51). Thus the time optimal synthesis for Σ S = [1, 2] is a combination of the individual synthesis for Neurons 1 and 2.

Spike Sequence [2, 2]
Similar controllability properties also allow us to give a short solution for the sequence Σ S = [2,2]. The optimal synthesis is shown in Fig. 4(c). In this case, the terminal cost ϕ(v 1 (τ 1 )) is a function of v 1 , and it is also strictly increasing in v 1 (since the higher the value of v 1 , the higher the time to ensure selective spiking in Neuron 2). From the analysis of transversality condition and the switching function like in the previous sequence (54) we can show that OFF-BANG is optimal for the first spike in Neuron 2 with a 1 < a 2 and suboptimal for a 2 < a 1 if there exists a switching. Indeed, for the first Neuron 2 spike and initial conditions under the separatrix, the optimal control is OFF-BANG. But for initial conditions on the v 2 axis, the optimal control is simply BANG. In the example, the overall construction is achieved by concatenating the solutions of P1 for Neuron 2 vertically. Since Neuron 2 is reset to 0 after firing, the initial condition for the second problem is given by

Proposition 3 Under Assumption 1, no time optimal control solution exists in general for a target sequence Σ S containing the subsequence [2, 1].
Proof The synthesis is more involved for this sequence. The terminal cost for the first Neuron 2 spike is the value function from (30) with v 2 = 0, that is, which is a decreasing function in v 1 , and ϕ(v 1 (τ 1 )) is not differentiable with respect to v 1 for some v 1 = v nd where v nd ∈ [0, V G ] (as shown in the bottom left of Fig. 5).
Note that for any initial condition at the origin or on the v 1 axis to the left of the separatrix, OFF-BANG cannot lead to optimality, and for those cases, the extremals will be generated by u * (t) = U for all t ∈ [0, τ 1 ]. Also, to the right of the separatrix OFF-BANG will be the optimal policy as it is the only viable option in the presence of state constraints. So we can conclude that if there is indeed a switching to the left of the separatrix, then there must exist v s ∈ (0, the optimal control is BANG. Now we will calculate this voltage v s , which acts as an onset for the change in optimal policy. Considering the switching at t =t, we have v 2 (t) = v s and Since the Hamiltonian vanishes identically for our problem, we get Also, from the transversality condition with λ 0 = 1 we have which is known. Since we reach the threshold V T from v s using the BANG control, from (26) we have Using the fact that the adjoint variables satisfy linear homogeneous differential equations, we can write From (59)-(63) we can solve for v s with If such v s exists, then the construction may be much more complicated with the possible presence of a "cut-locus" type phenomenon, and we leave a detailed analysis of such a problem for future work. In our case, the terminal cost decreases with a rapid rate for v 1 ∈ [0, v nd ] and abruptly changes to a much smaller slope for v 1 ∈ (v nd , V G ] (see Fig. 5) due to the nature of the value functions on either side of separatrix V−, V+ in (31) and (32). This results in a field of extremals trying to converge to the point v nd , even when the monotonicity of the value function is not affected by the loss of differentiability (see top left in Fig. 5). We calculate the set of initial conditions for which this point can be attained, where v c denotes the highest point on v 2 axis from which [V T v nd ] T can be reached via BANG control. This voltage v c and the set v c are shown in the right panel of Fig. 5. Now, the optimal control problem for v(0) ∈ v c simply reduces to This is similar to the selective spiking problem of Neuron 1, and indeed the best control is a combination of BANG and boundary control as in (17), But this implies that Neuron 2 maintains the voltage (V T ), even after the spike is emitted, which violates our assumption that the neurons are reset instantaneously after reaching V T , as described in (4). So the synthesis S * corresponding to (66) is excluded from the admissible set of extremals purely out of the physical constraints imposed on the system. This resembles the classical problem of finding surfaces of minimum revolution [34], where the Goldschmidt extremal cannot be attained because of the C 1 assumption on the extremals. Thus, any synthesis S for (65) will be suboptimal to S * . For simplicity, we have picked a synthesis such that that is, a constant control that varies depending on the initial condition shown in Fig. 5. For the set of initial conditions the optimal synthesis remains the same as the solution of P1 for Neuron 2.

Greedy Designs for Sequences with Arbitrary Length
From our analysis of the 2-spike sequences in the previous section, we can design the time optimal control for any Σ S of M spikes (M ≥ 2) without the subsequence [2,1]. In addition, if we assume that a 2 < a 1 , then it can be shown using an inductive argument that the overall synthesis can be constructed from the solutions of individual selective spiking problems in Propositions 1 and 2.
In general, for a Σ S with the subsequence [2,1], to illustrate the complexities of sequence control, it is instructive to consider the 4-spike sequence Σ S = [1, 2, 1, 1]. In this case, the target sequence contains a [2, 1] event, meaning that any solution will be suboptimal. In this case, a dynamic programming approach that interleaves the interpolation control (67) can yield such a solution. However, from a practical perspective, pursuing this design approach for long sequences is difficult as it requires computing the location of nondifferentiability in the value functions of all [2, 1] events. Fig. 6 Simulation example of the greedy algorithm for P2 for a target sequence Σ S = [1, 2, 2, 1, 1, 2, 1] with the nominal parameters in (39). The inset shows the synaptic contribution v 2 (t) = 2 mV, to Neuron 2 due to the first spike in Neuron 1.
Thus, we argue that, from a design perspective, a simple greedy approach, where we minimize the time for each spike in Σ S progressively, constitutes an acceptable and tractable approximation.
In Fig. 6, we show the solution of the greedy controller for an arbitrary spike sequence Σ S .

Decoupling the Network for Longer Sequences
In applying the greedy approach, it is important to note that the synaptic contribution from the spiking neuron can carry the voltage of the other neuron in the network over the synaptic guard V G . Thus, we cannot readily apply the solution of P1 for the following spike in the sequence (pattern), as the initial condition may violate the state constraint in (8) for P1. Here, we propose strategies to eventually utilize Propositions 1 and 2 for the greedy design.
1. First, if the initial condition after any spike in the sequence (pattern), at t = τ 1 , is not within the relevant state space G, then we can apply u = 0 until t = t , t > τ 1 , such that v(t ) ∈ G. Then, we can apply the solution of P1 to induce the target spike. 2. Alternatively, we can modify the guard V G of the nontarget neuron at each step of the greedy design, depending on the number of consecutive spikes in the target neuron in the sequence (pattern); for example, if Σ S = [1, 1, 2, 2, 2], then we can set the guard voltage for Neuron 2 at V G (σ 1 ) < V T − 2ρ syn for the first spike and V G (σ 2 ) < V T −ρ syn for the second spike. Thus, the relevant state space for the first and second spikes will be modified to G( respectively. This ensures that whatever the contribution is from the presynaptic neuron (in this case, Neuron 1), we start in the relevant state space for the next spike in the sequence (pattern). Once the target neuron changes to σ 3 = 2, the guard voltage for Neuron 1 is determined by the number of consecutive spikes in Neuron 2 (3 in this example), that is, V G (σ 3 ) < V T − 3ρ syn and so on. Note that by successively reducing the guard voltage, the selective spiking problem may become infeasible as discussed in Sect. 3.3. 3. Finally, we can combine the two approaches to develop an algorithm where we can use (2) until the problem is infeasible. At this point, we go back to (1) and add an off time before implementing the solution of P1.
In our examples of sequence and pattern control, we have used the first approach in developing the greedy design (see Figs. 6 and 7).

Fixed-Time Selective Spiking and Spike Patterns
We now move to the problem of controlling timed spike patterns, that is, P3. It is intuitive that a basic necessary condition in this case is that the desired spike time exceeds the minimum selective spiking time, that is, the solution to P1. Specifically, suppose that we want to achieve the target pattern Σ P = [(1, t 1 )], that is, a spike in Neuron 1 at time t 1 . The cost function in P3 (14) reduces to (subject to the selectivity constraint in (7)). Here, τ 1 denotes the achieved spike time, andτ 1 is the solution of P1 for an arbitrary initial condition v(0). Ifτ 1 ≥ t 1 , then evidently that is our best option, and the solutions of (69) and P1 are the same, that is, τ 1 =τ 1 .
For the other case,τ 1 < t 1 , contingent on controllability, a control must exist such that τ 1 = t 1 . If such a condition is met, then in general there may be multiple solutions to the pattern control problem.
Herein, we consider one simple strategy involving the introduction of an off timê t to the optimal control solution of P1 such that where τ r 1 is the solution of the time optimal control P1, for the initial condition v(t). We noted earlier that the initial conditions for the selective spiking problem nominally lie on either the v 1 or v 2 axis, under the assumption that one of the neurons has just produced a spike. In this case, feasibility of (70) reduces to understanding those initial conditions that generate specific values of τ r 1 .

Off-Time Insertion for Pattern Control
We characterize the relationship between τ r 1 and initial conditions via the notion of a Λ-controllable set.
Definition 3 (Λ-Controllable set) Without loss of generality, the Λ-controllable set ζ(Λ) of Neuron 1 is the set of initial conditions from which the selective spiking of Neuron 1 in P1 is achieved in time Λ, that is, The Λ-controllable sets for system (11) are provided in Appendix A.4. Since we are interested in initial conditions along the v 1 and v 2 axes, we consider the functions that is, the intersection of the Λ-controllable sets with the axes.
Earlier, we noted that the value function for the selective spiking of both neurons remains continuous on both the v 1 and v 2 axes (i.e., from (36) and (44)). This fact, together with the derivation of the Λ-controllable sets in the Appendix, allows us to conclude that the functions (72) are continuous in Λ.
Thus, we are able to ensure the existence of the off-time pattern control from (70), that is, where u * comes from Proposition 1 or 2. The computation of the off-timet is obtained directly from the Λ-controllable sets and is provided in Appendix A.5. Thus, an overall pattern control strategy can be formulated as

Greedy Designs for Control of Long Patterns
We now consider the synthesis and design of the general pattern control problem P3.
To begin, we consider the dynamic programming strategy studied in (51) but for P3. It turns out that the same issues pertaining to nondifferentiability of the value function in P2 persist in this case.
To illustrate this, consider the 2-spike target pattern Σ P = [(1, t 1 ), (1, t 2 )]. Starting from the last spike σ 2 = 1, we solve with the terminal and state constraints and use the value function of (75) as the terminal cost to the following optimal control problem: Let us denote the solution of P1 for the second spike from the initial condition v(0) = [0 v 2 ] T byτ . Then, depending on v 2 , the terminal cost in (76) takes the following form: Thus, similar complications as referenced in Sect. 4.2 regarding nondifferentiability arise here, and once again we consider implementation of a straightforward greedy strategy for pattern control involving (74). In Fig. 7, we show an example of this greedy algorithm for an arbitrary pattern with the same spike sequence as in Fig. 6.  (39). Similar to Fig. 6, we show the synaptic contribution v 2 (t) = 2 mV, to Neuron 2. We also explicitly indicate the off-time (u = 0) after the first (inset) and fourth spike in Neuron 1, as part of the decoupling strategy discussed in Sect. 4.2.1.

Performance of Greedy Design Under Disturbance and Noise
In this section, we analyze the robustness of the greedy design when the coupled LIF network in (3) is subjected to noise and disturbances. Here we consider two types of uncertainties: 1. Incoming synaptic contributions of the pulse coupled form discussed in Sect. 2.2.2, from other neurons 2. Noise in the dynamics of the membrane voltage of the neurons in (3) (process noise) and in measurement of these voltages (measurement noise). Note that in implementing the greedy controller in (74), we repeatedly apply Propositions 1 and 2, which are feedback control, that is, measurement is implicit.
In Fig. 8(A), we show one realization of the voltage and control waveforms for d = 150 incoming spikes over the control horizon for the same Σ P used in the example of Fig. 6. To illustrate the effect of these disturbances on the control strategy, in Fig. 8(D), we plot the average Victor-Purpura (VP) distance [36,37] between the achieved and target spike trains as we vary the number of incoming spikes d over 50 trials. In each trial, we randomly select the arrival times of the spikes, the contribution and target of the synapse between the two neuron indices. The VP metric is a measure of synchrony between two spike patterns that involves three basic operations: adding or deleting any spike with cost 1, moving any spike with cost q per unit time, and renaming any index of the spike with cost k. Here, a lower VP distance corresponds to better control performance. We observe that with higher disturbance, represented by d, the controller performs reasonably well with gradual degradation in the achieved patterns. Next, we consider additive Gaussian noise both during the evolution of the membrane voltage and in measurement. Thus the linear model in (11) is modified tȯ v(t) = Av(t) + bu(t) + w(t), where the measurement vector y is a linear readout of the neuron voltages through a randomly selected matrix C, which is full rank, w(t) and z(t) follow multivariate Gaussian distributions with w(t) ∼ N (0, W ) and v(t) ∼ N (0, Z), and W and Z are the constant covariance matrices of the forms W = η 2 1 I and Z = η 2 2 I, I is the 2-dimensional identity matrix. Here, we compute the voltage estimates of the two neurons at each time step by means of a Kalman filter [38] and employ the feedback strategy in (74) based on these estimates. In Fig. 8 (B), (C), we plot the pattern control solutions for the same Σ P used in the example of Fig. 6 for smaller (η 1 = 0.1, η 2 = 1) and higher (η 1 = 1, η 2 = 10) process and measurement variance. We observe that controller's ability to induce the target spike train is not compromised substantially, although with higher levels of noise, spurious spikes are generated, as indicated in panel (C). However, the noisy dynamics in (78) can result in a high frequency of switching in the control ((B), (C), bottom panel), especially during the boundary arc, that is, the nontarget neuron is to be held at guard V G . Panel (E) shows the performance of the greedy design with respect to the average VP distance between Σ P and achieved patterns over 50 different trials, as we change the level of noise during the evolution and measurement phase.

Conclusions
This paper has examined the problem of controlling timed spike patterns in pairs of Integrate and Fire neurons. Boundary-arc-type phenomena are shown to arise in this scenario due to state constraints imposed by both the selectivity criterion and spike generation mechanism. Formal analysis and synthesis is carried out to establish how the proposed solutions are geometrically disassociated in terms of their initial conditions. The developed solutions, which leverage the maximum principle and dynamic programming, are shown to be efficacious in controlling the LIF models.
Clearly, our results here are of theoretical nature. Although the control-theoretic features revealed are themselves interesting from a mathematical standpoint, they serve the broader purpose of establishing fundamental limits on the selective control of neurons with common inputs. The qualitative nature of the derived solutions (e.g., OFF-BANG, boundary-arc strategies) are already more complex than the fixedamplitude square pulse designs currently used in practice. Given the massive growth in stimulation technology development, understanding these limits, even for a relatively simple model class, may provide insight into how experimentalists should tune their stimulation parameters for experimental objectives. For instance, our analytical conditions (e.g., γ 1 ≷ V T V G ) amount to a criteria on the amount of heterogeneity needed within a neuronal population in order to enable control. Without sufficient heterogeneity, it is simply impossible for a common input to "split" the spiking of neurons in a selective manner. Exploiting this heterogeneity is at the heart of the derived control solution (e.g., OFF-BANG solutions that leverage increased leak dynamics). These characterizations provide a baseline from which we plan to establish relaxations of the considered problems for larger neuronal populations in future work. with where t s is the arrival time of a presynaptic action potential from the other neuron, g syn is the synaptic conductance, H is a Heaviside step function, κ s is the time constant for the conductance,ḡ s is the maximum conductance for the synapse, and E syn denotes the reversal potential. For selective spiking, we want the postsynaptic neuron to be protected from this incoming synapse with respect to the membrane potential.
In the typical case of an excitatory synapse, we have E syn ≈ 0, and the contribution from the spike in the presynaptic neuron becomes Now, assuming a separation in time scale between the synaptic time constant κ s and the membrane time constant κ m , that is, κ s κ m , we can approximate the integral in (81) by keeping the voltage of the postsynaptic neuron constant at v(t s ) over the integration window. Using this, we have So the effect of a synaptic event on the postsynaptic neuron can be crudely summarized as an almost instantaneous rise in voltage bounded by (82). Thus, model (5) approximates this effect with an impulsive synaptic action, where

A.2 Geometrical Aspects of Selective Spiking Solution
Here, we first discuss the role of γ 1 (15) in determining the two different selective spiking solutions presented in Sects. 3.1 and 3.2. We also geometrically show that pairwise feasibility is not achievable when both neurons are Case 2, as described in Sect. 3.3.2. We first derive the equation for the line of quasistatic equilibrium defined in (48). This is the set of points in the phase plane for whichv(u) = 0 for any constant control u ∈ U . Using this condition, we havė Since u is a constant, we can eliminate u to get the equation for the quasistatic equi- where γ 1 = b 1 a 2 b 2 a 1 . Now the two different solutions presented in Propositions 1 and 2 are dependent on the existence of the boundary segment, that is, for a boundary control u arc for which Neuron 2 is voltage invariant (v 2 (u arc ) = 0), regardless of whether the voltage of Neuron 1 increases. To satisfy this, we must havė We can answer this question from the analysis on the quasistatic equilibrium line as u arc is constant. If the line intersects v 1 = V T before v 2 = V G , from (85) we have (see Fig. 9(a)) Using this, we can calculate the direction of vector field at v 1 = V T in (86): We show this in Fig. 9(a), where the quasistatic equilibrium intersects v 2 = V G beyond v 1 = V T . This ensures that the vector field is positive under the boundary control such that the target neuron reaches threshold while keeping the other neuron at V G . Now if we assume that γ 1 ≤ V T V G (see Fig. 9(b), (c), (d)), then we can similarly show as in (88) thatv for this case, and we need to adopt the solution presented in Proposition 2 to fire Neuron 1 selectively. So we see that the nature of selective spiking solution, that is, (BANG/BANG-BOUNDARY) or (BANG/OFF-BANG), is contingent upon the ratio γ 1 . Figure 9 presents an intuitive representation of the geometric aspects of the solution space discussed in Sects. 3.1-3.3 with respect to γ 1 , γ 2 . Here, we analyze the pairwise feasibility for all possible parameter combinations. If γ 1 > V T V G (=⇒ γ 2 < V T V G ) and Lemma 1 for Neuron 2 holds, then the neurons are pairwise feasible, that is, from any point in the phase plane we can fire either neuron selectively. Similarly, if we have , that is, Neuron 2 is Case 1, Neuron 1 is Case 2, and Lemma 1 holds for Neuron 1, then we can once again achieve pairwise feasibility. These two scenarios are depicted in Fig. 9(a), (d), respectively.
, both neurons are Case 2), for pairwise feasibility, we must have Lemma 1 satisfy for each neuron individually. This creates a situation shown in Fig. 9(b), (c), where the separatrices for Neuron 1 and Neuron 2 intersect, which implies that, at the point of intersection, we have two different vector fields under the same control (u = U ), which is a contradiction. Hence, if both neurons are Case 2, then we cannot have pairwise feasibility.

A.3 Sketch of Regular Synthesis-Type Sufficient Conditions for Optimality
We outline regular synthesis-type sufficient conditions for optimality and how they apply to the optimal control problems considered in this paper. Any sufficiency theory in optimal control problems is based on the fact that the cost-to-go function V of the optimal control problem evaluated along extremal controls u * in a synthesis (to be verified as the optimal one), that is, the value of the objective is evaluated along the controls u * as a function of variable initial time t 0 and initial condition v 0 , is, in a suitable way, a solution to the Hamilton-Jacobi-Bellman equation (38), If the cost-to-go function is continuously differentiable everywhere-and this is Bellman's classical argument-then this is easy: given any other controlũ with corresponding trajectoryṽ defined over an interval [t 0 , T ], it follows that and thus Hence and thus the cost J (ũ) alongũ is not better than the cost along the control u * given by V(t 0 ,ṽ 0 ). This proves the optimality of the controls in the synthesis.
The main issue with this argument is that problems for which the cost-to-go function is differentiable everywhere are generally rare. It is a much more common scenario that the value function loses differentiability on thin subsets, typically given by locally finite unions of embedded submanifolds of positive codimension. For the problems considered in this paper, such examples are given by the separatrix in Example 2 or a cut-locus that arises for the spike sequence [2,1]. In fact, such structures are omnipresent in solutions of time-optimal control problems.
This argument breaks down if there exist lower-dimensional submanifolds along which V is not differentiable since, in principle, the set of times when the controlled comparison trajectoryṽ lies in such a submanifold can be an arbitrary closed subset of the interval [t 0 , T ], and it is simply no longer possible to differentiate the function V along such a trajectory. Alleviating this issue is a highly nontrivial technical matter, which has led to regular synthesis-type arguments for optimality [39].
The key technical step is to perturb the given comparison trajectoryṽ in such a way that the perturbed trajectory has a cost that is close to J (ũ), whereas it meets the manifolds where V is not differentiable only in a finite set of times. It is clear that then the previous argument can be carried out piecewise, and the result follows in the limit as the approximations approachṽ. We briefly sketch the main steps of this reasoning, but refer the reader to Chap. 6.3 in [34], where the technical details of this argument are carried out in full.
In the first step, the comparison controlũ (a priori only a Lebesgue-measurable control) is approximated by a sequenceũ n of piecewise constant controls for which both the corresponding trajectoriesṽ n converge toṽ uniformly over the interval [t 0 , T ] and the integrals converge as well, Note that the approximating controls are piecewise constant, that is, have only a finite number of switchings, and are not just simple measurable controls with a finite number of values as in the definition of Lebesgue measurability.
In the second step, the initial time t 0 and initial condition v 0 are perturbed so that the resulting trajectoriesṽ n only meet the submanifolds where the cost-to-go function V is not differentiable in at most a finite number of points. Essentially, this is just a transversality argument: since the perturbed controls are piecewise constant, trajectories are integral curves of smooth vector fields. If the flow of a vector field at a point p is transversal to an embedded submanifold M with positive codimension, then there exists ε > 0 such that the flow does not lie in M for 0 < |t| ≤ ε, and thus the corresponding time is an isolated point of the set of times when the flow meets M. We just need to make sure that there exist enough initial conditions for which this holds, and this can be guaranteed using Sard's theorem [40].
As a consequence of all these perturbations, we need to take the limit as n → ∞ and the trajectories generally no longer satisfy the terminal constraint. But all this works out in the limit if the cost-to-go function V is lower semicontinuous (and this is a necessary condition for it to be the value function of an optimal control problem in which the objective is minimized) and in addition satisfies the following two technical requirements: (1) For every constant admissible control u, the function V has the no-downwardjumps property along the vector field x → f (x, u), that is, if γ is an integral curve of such a vector field defined on a compact interval [a, b], γ : [a, b] → G, t → γ (t), then, for all s ∈ (a, b], we have that (2) For every point q in the terminal manifold and every ε > 0, there exists a nonempty open set Ω ⊂ G ∩ B ε (q) such that, for all z ∈ Ω, we have that These two conditions are satisfied automatically wherever V is continuous.
For the problems considered in this paper, in Case 2 the cost-to-go function V need not be continuous, but it is not difficult to argue that condition (1) is satisfied. The value cannot decrease along constant controls as the separatrix would be crossed. This simply holds since trajectories can only cross from the region where the cost is smaller into the region where the cost is higher as the separatrix is defined by the maximum value of the control. Similarly, condition (2) is a weak continuity requirement that allows discontinuities in the value function at the terminal manifold. It only requires that, for any potential target point q, there exist sufficiently rich sets that are close such that the cost-to-go function V still has some upper continuity property along sequences converging to q in those sets. For our problems, this is satisfied by simply choosing Ω to lie in the regions where the optimal control is at its maximum value. Thus these conditions, albeit technical, are quite natural and generally are easily verified as is the case for our problems.
We close with some remarks about the synthesis in Case 1 when the trajectories contain the boundary arc. In this case, following the guard is in fact the only feasible control that the system can take, and thus, once the system hits the guard, the only control possible is the boundary control. This allows us to eliminate the state space constraint by extending the terminal set to include the guard v 2 = V G . We only need to define a penalty term ϕ 2 on the guard that represents the time it takes for the boundary control to reach the terminal state v 1 = V T . This transforms the time optimal control problem with state space constraint v 2 ≤ V G into an unconstrained optimal control problem whose terminal set consists of the union of two manifolds: the regular terminal manifold v 1 = V T with penalty ϕ 1 ≡ 0 and the guard v 2 = V G with penalty function ϕ 2 . The two manifolds intersect at the point [v 1 v 2 ] T = [V T V G ] T , and ϕ 2 = 0 at this point. In this case, the cost-to-go function is not differentiable along the trajectory that ends in this point, but the conditions mentioned before are satisfied, and this is the optimal synthesis.
We emphasize that this argument is merely a convenient trick, which becomes available for this particular situation in Case 1. In general, however, that is, for Case 2 and the subsequent time optimal control problems of spike sequences, the methods and techniques sketched before are essential. Standard sufficiency arguments based on viscosity solutions to the Hamilton-Jacobi-Bellman equation fail once the costto-go function becomes discontinuous. This feature, related to questions of smalltime local controllability (e.g., see [34,Sect. 7.2]), is no obstacle for synthesis type sufficiency arguments.

A.4 Computation of Λ Controllable Sets
We now show the calculation for Neuron 1. There are two possible situations, namely, Λ ≤ T s and Λ > T s , which result in two different switching structures where T s denotes the time to reach [V T V G ] T along the separatrix Γ from the initial condition If Λ ≤ T s , then we can find the neuron voltages (v 1 , v 2 ) from which Neuron 1 reaches V T , in time Λ, Note that v 2 does not come in (98) since for all v ∈ Γ + , Neuron 1 reaches threshold without Neuron 2 hitting the guard. For Λ > T s , we assume that it takest for Neuron 2 to hit the guard V G under bang control, The voltage of Neuron 1 at this time is calculated using (33). This means for (v 1 , v 2 ) to be on the Λ-controllable set, Neuron 1 must reach the threshold V T in (Λ −t) along the boundary arc, that is, Simplifying (100), we get where a 2 . From this we can find the Λ controllable set for the selective spiking of Neuron 1. Similarly, for Neuron 2, we can find the set ζ 2 (Λ).

A.5 Calculation of off-Time for Fixed-Time Selective Spiking
In this section, we show how the off-time in (70) can be calculated to induce a spike in a specified time. Without loss of generality, we once again assume the target pattern . (102) , we get the desired off-time. Note that, for t 1 T s , we will need to use the boundary segment in (102).

Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.