Ill-Posed Point Neuron Models

We show that point-neuron models with a Heaviside firing rate function can be ill posed. More specifically, the initial-condition-to-solution map might become discontinuous in finite time. Consequently, if finite precision arithmetic is used, then it is virtually impossible to guarantee the accurate numerical solution of such models. If a smooth firing rate function is employed, then standard ODE theory implies that point-neuron models are well posed. Nevertheless, in the steep firing rate regime, the problem may become close to ill posed, and the error amplification, in finite time, can be very large. This observation is illuminated by numerical experiments. We conclude that, if a steep firing rate function is employed, then minor round-off errors can have a devastating effect on simulations, unless proper error-control schemes are used.


Introduction
Modeling of electrical potentials has a long tradition in computational neuroscience. One model with some physiological significance is the voltage-based system τ u (t) = −u(t) + ωS β u(t) − u θ + q(t), t ∈ (0, T ], where u(t), q(t) ∈ R N , t ∈ (0, T ], In the rate model (1)- (2), each component function u i (t) of u(t) represents the time dependent potential of the ith unit in a network of N units. The nonlinear function S β is called the firing rate function, {ω ij } are the connectivities, and q(t) models the external drive. A detailed derivation of this model can be found in [1][2][3]. The purpose of this paper is to explore the properties of the initial-condition-tosolution map associated with (1)- (2). Note that we use the subscript β to emphasize that R β depends on the steepness parameter β, and that R ∞ corresponds to using a Heaviside firing rate function, i.e. S ∞ = H . We will also make use of the standard notation for the supremum norm throughout this paper. A simple example, presented in Sect. 4, shows that R ∞ can become discontinuous. Hence, the model is mathematically ill posed [4,5] and round-off errors of any size can corrupt computations. We conclude that it is very difficult to produce reliable simulations with such models. Since all norms for finite dimensional spaces are equivalent, it is not possible to "circumvent" this problem by changing the involved topologies.
According to standard ODE theory (Appendix A), R β , with β < ∞, is continuous, but the size of the error-amplification ratio may be huge for large β, which will be demonstrated and analyzed in Sects. 2 and 3, respectively. Here,ũ 0 represents a perturbed initial condition andũ(t) its associated solution. This implies that, also for 1 β < ∞, it can become difficult to guarantee the accurate numerical solution of (1)-(2): Minor round-off errors may be significantly amplified within short time intervals, which can lead to erroneous simulations. Our investigation is motivated by the fact that steep sigmoid functions, or even the Heaviside function, often are employed in mathematical/computational neuroscience; see e.g. [1,6] and references therein. Other authors [7,8] have also pointed out that severe challenges occur if β = ∞, i.e. issues concerning how to define suitable function spaces and to prove existence of solutions. Nevertheless, as far as we know, results which explicitly discuss the ill-posed nature of (1)-(2) when β = ∞, and how this property yields extra numerical challenges in the steep, but smooth, firing rate regime, has not previously been published.
Remark We would like to point-out the following: Assume that an initial condition is close to an unstable equilibrium. Our results should not be interpreted as expressing the mundane fact that a perturbation of this initial condition, moving it to another region with completely different dynamical properties, may lead to large changes in the solution. In fact, we show that the error-amplification ratio can be huge, during small time intervals, even though the perturbation does not change which neurons are active. That is, the change in the initial condition is not such that it changes the qualitative behavior of the dynamical system for 0 < t 1-only the quantitative properties are dramatically altered. This can happen in the steep firing rate regime.

Numerical Results
Let us first compute the error-amplification ratio (5) for some simple problems.
Example 1 Consider the following model of a single point neuron, i.e. N = 1, We used Matlab's ode45 solver, with the default error-control settings and T = 0.1, to compute numerical approximations of Also a second series of simulations were performed, using the same selection of values for the steepness parameter, but with the perturbed initial conditioñ The corresponding solution is denotedũ(t) =ũ(t; β). Plots of u andũ, with β = 200, are displayed in Fig. 1. Note that in both cases the neuron fires, i.e. the change in the initial condition is not such that it has moved from one side of an unstable equilibrium to the other side. Even so, according to Fig. 2 and Table 1, the error-amplification ratio E(T ; β), due to the minor perturbation (6) of the initial condition, is in the range [80.6, 1054.1] for β ∈ [100, 200], and also very large for β = 50, 75.
Simulations with the strict error-control setting odeset('RelTol',1e-13,'AbsTol',1e-13) generated the same results, and so did an explicit Euler scheme, with uniform timestep t = 10 −7 .  Example 2 Let us consider a model of two point neurons:  The same procedure as in Example 1 was used, but with the perturbed initial conditioñ produced virtually the same results. The simulations were also "confirmed" by our explicit Euler implementation with time-step t = 10 −7 . Figure 6 shows numerical results computed with Matlab's ode15s solver, employing the default error-control settings. The curves shown in this figure are very different from the graphs displayed in Fig. 4, which were computed by the ode45 software. We conclude that even the toy example considered in this section is not  trivial to solve (with the strict error-control setting (7), ode15s also managed to produce the curves shown in Fig. 4). If and the model implies that which are rather small. One therefore might think that it is sufficient to employ a moderate time-step to obtain an accurate numerical approximation. Figure 7 shows that this is not the case. (In computational mathematics it is well known that the accuracy of the finite difference approximation [u 1 (t + t) − u 1 (t)]/ t, of the derivative u 1 (t), depends on the second order derivative u 1 (t), which in our case is of order O(β). This explains the poor approximation obtained with time-step t = 0.01.)

Analysis
The purpose of this section is to present an analysis of the error-amplification ratio (5) and thereby explain the main features of our numerical results. Even though the Picard-Lindelöf theorem [9,10] asserts that (1)-(2) has a unique solution u(t), provided that q(t) is continuous and that β < ∞, it is virtually impossible to determine a simple expression for u(t). On the other hand, if u(t) ≈ u θ and β < ∞, then we can linearize S β to get an approximate model, which is much easier to work with.

Linearization
The linearization of S β about zero reads Define τ = I, the identity matrix, then the linear approximation of (1)-(2) reads where The linearized problem with a perturbed initial condition becomes Therefore, and the error-amplification ratio can be written in the form Since the entries of A = A(β) are of order O(β), see (11), we conclude that the erroramplification ratio for the linearized model is of exponential order O(e β ). Is this also the case for the highly nonlinear model (1)-(2)? We will now explore this issue, but first we would like to make a short remark.
Remark Recall the definition (8) then the analysis of the linearized model, presented above, would also be valid for (1)-(2), provided that Similarly to the sigmoid function S β ,L β also converges point-wise to the Heaviside function as β → ∞. If one employs the sigmoid function in the point-neuron model, then the analysis, as we will see below, becomes much more involved.

Preparations
Let β max ,T , and α be arbitrary positive constants. It is easy to construct a smooth vector-valued function z satisfying Hence, defining the source as we conclude that the solution u(t; β max ) = z(t) of (1)-(2) also satisfies (13), provided that u 0 = z(0). By employing standard techniques, one can show that the solution u(t; β) of (1)-(2) depends continuously on 0 < β < ∞; see Appendix B. Consequently, there existsβ min < β max such that For the sake of simple notation, we will in our analysis write u, or u(t), instead of u(t; β). Furthermore, according to the analysis presented in Appendices A-C, u depends continuously on both the initial condition u 0 and the steepness parameter β, when 0 < β < ∞. Motivated by this property of (1)-(2), we assume that both u andũ, whereũ denotes the solution of (1) generated by a perturbed initial conditionũ(0) = u 0 , satisfy whereβ min < β max . Then, by invoking the triangle inequality, we find that which will be small if β max is large. Even so, as will become evident below, the error-amplification ratio (5) can be significant and lead to erroneous results. Let s ands denote the associated solutions of the linearized model (9)-(10). From (14) we find that the initial conditions u 0 andũ 0 satisfy Since s ands are continuous with respect to t, the same initial conditions are employed in the linearized model, and these solutions depend continuously on 0 < β < ∞, it follows that there existT > 0 andβ min < β max such that The main point of this discussion is to show that there exist (smooth) source terms q and perturbations of the initial condition such that (14) holds, regardless how largê T , β max , α > 0 are. Also, the solutions of the linearized model will satisfy (16). For the sake of simple notation, let T = min{T ,T } and β min = max{β min ,β min }.
The triangle inequality implies that We will derive a bound for e(T ) ∞ . The analysis of ẽ(T ) ∞ is completely analogous, and thus it is omitted.

Linearization Error
Subtracting (9) from (1), and keeping in mind that we consider the case τ = I, yields The triangle inequality, Taylor's theorem and Eq. (8) for L β imply that where the second last inequality follows from (14). By combining this with (18), and the triangle inequality, one finds that Since this must hold for i = 1, 2, . . . , N, and Grönwall's inequality implies that

Error-Amplification Ratio
Clearly, and the reverse triangle inequality yields From (12) it follows that the error-amplification ratio (5) satisfies
Recall that the entries of the matrix A = A(β) are of order β; see (11). To derive a bound for II(T ; β), we employ (19), and a similar inequality for ẽ(T ) ∞ , is not very large, β is fairly large and, e.g., α ≥ 0.5, then the size of the erroramplification ratio E(T ;   Remark Assume that the · ∞ -norm of the source term q(t) is bounded. Then, since 0 < S β [x] < 1 for all x ∈ R, it follows from (1) that both u (t) ∞ and ũ (t) ∞ are bounded independently of the size of the steepness parameter β, at least when u(t) ≈ u θ andũ(t) ≈ u θ . Consequently, also the difference u(T ) −ũ(T ) ∞ is bounded independently of β > 0. Our results therefore might appear to be somewhat counterintuitive: But note that we have only argued that the error-amplification ratio (5) may, approximately, be of order O(e β ). If β is large, this can cause severe numerical challenges.
We would also like to comment that standard theory for general dynamical systems relies on the size of F , which for the point-neuron model (1) is not Lipschitz continuous with respect to z when β = ∞, which the Picard-Lindelöf theorem [9,10] requires. (F is not even continuous when β = ∞.) The maximum error bound (15), valid when u − u θ andũ − u θ satisfy (14), suggests that setting β = ∞ might provide a solution to the issues discussed above. Unfortunately, as will be explained in the next section, this is not the case.

Ill Posed
We will now show that (1)-(2) can become truly ill posed, if a Heaviside firing rate function is employed. More specifically, the initial-condition-to-solution map, in finite time, can be discontinuous.
Consider the case N = 1, τ = 1 and no source term: If, for 0 < 1, provided that ω > u θ . Consequently, where R ∞ denotes the initial-condition-to-solution map (3). We conclude that, no matter how close u 0 and u 0 are, the difference v(T ) − v(T ) between the corresponding solutions will not become small. Hence, R ∞ is discontinuous. It follows that the initial value problem, with a Heaviside firing rate function, is ill posed, in finite time-at least in the sense of Hadamard. Also note that, unless ω = 2u θ , u θ is not a stationary solution of (22), i.e. not an unstable equilibrium. The error-amplification ratio for this ill-posed problem becomes infinite when → 0: as → 0, for any T > 0. One may consider this issue from a more pragmatic point of view. Let v t denote a numerical approximation of v. If a Heaviside firing rate function is employed, then H (v t − u θ ) must be evaluated in some line of the simulation software. This is an unstable procedure because H has a jump discontinuity at 0, and round-off errors of any size can corrupt computations.
In contrast to this, provided that β < ∞, see the analysis of the model (1)-(2) presented in Appendix A. Here,ũ 0 is any perturbation of the initial condition u 0 , and A and B are positive constants depending on the matrices τ and ω, but not on β. This inequality shows that the initial-conditionto-solution map R β , β < ∞, also is continuous at unstable equilibria.

Conclusions and Discussion
Since R ∞ can become discontinuous, it is virtually impossible to guarantee the accurate numerical solution of point-neuron models which employ a Heaviside firing rate function: Any round-off errors can potentially corrupt simulations. Alternatively, one may stop the simulation as soon as the solution hits the jump discontinuity, i.e. the threshold value for firing. We have also observed that models with a steep, but smooth, firing rate function can amplify errors to an extreme degree, which is typical for "almost ill-posed" problems. Consequently, reliable simulations can only be obtained if proper error-control schemes are invoked. How to design effective error-control methods, for models with a large steepness parameter β, is, as far as the authors know, still an open problem. Nevertheless, it seems plausible that suitable adaptive numerical schemes, where the time steps become smaller when the solution reaches regions in the vicinity of the threshold value for firing, might be capable of handling the numerical error amplification. Let be the operator which maps the solution of the point-neuron model (1)-(2) from time t 1 to time t 2 . Note that the action of F β;t 1 ,t 2 can be determined by solving the pointneuron model with u(t 1 ) as initial condition. Therefore, from the argument presented above, it follows that the error amplification ratio associated with F β;t 1 ,t 2 may be large, provided that β 1. We conclude that the issues pointed out in this study cannot necessarily be avoided by using an initial condition which is far from the threshold value u θ for firing. In fact, it seems that one must prove that u(t) never gets close to u θ for t > 0-a herculean task, if correct.
From a modeling perspective one might wonder: Should a voltage-based model of cortex be ill posed or "almost ill posed"? If so, then models employing a Heaviside firing rate function cannot be robustly solved with finite precision arithmetic and regularized approximations are numerically challenging [4,5].
We fear that similar unfortunate properties, to those discussed in this paper, might be valid for models which can be written in the form where F β → ∞ when β → ∞. This can, e.g., be the case for a number of models in use in computational neuroscience and gene regulatory networks.
An easy solution to the issues raised in this paper, is to avoid steep firing rate functions. If β is fairly small, then standard ODE theory [9,10] and textbook material about their numerical treatment can be used, provided that the source term q(t) is continuous. Nevertheless, steep sigmoid functions are popular in computational neuroscience.

Appendix A: Continuous Dependence on the Initial Condition
We will prove the stability estimate (23), which holds if β < ∞. Let u andũ denote the solutions of (1) corresponding to the initial conditions u 0 andũ 0 , respectively. For the purpose of detailing the derivation of this estimate, we work with the components of these vector functions, i.e., Let us also assume that the diagonal entries {τ i }, of the diagonal matrix τ , are positive, and let {ω ij } denote the entries of ω.
We first observe that the component functions u i andũ i satisfy the fixed point problems from which we find, by using the triangle inequality, that Then, by exploiting the fact that for all x, y ∈ R, we get where we in the final step have used the definition of the supremum norm (4). Thus, we arrive at the inequality The stability estimate (23), with then follows from Grönwall's lemma.

Appendix B: Continuous Dependence on the Steepness Parameter
For the sake of completeness, we also show that the solution of the initial value problem (1)-(2) depends continuously on the steepness parameter β. Let β andβ be steepness parameters for the firing rate function. We fix the initial condition and the connectivity parameters of (1)- (2). The solutions corresponding to β andβ are denoted by u andû, respectively. We readily obtain for the component functions u i andû i of u andû, respectively. We now make use of the property for the firing rate function, and we get the chain of inequalities |ω ij | |u j,θ |T + which proves that the solution of (1)-(2) depends continuously on the steepness parameter β < ∞ of the firing rate function.
Since q(T ) is increasing, it is possible to extract further information from this argument. More specifically, let U = v(t) v(t) ∈ R N , t ∈ [0, T ] and sup Then we can conclude that the mapping is continuous.

Appendix C: Continuous Dependence on the Initial Condition and the Steepness Parameter
We finally prove that the solution u = u(t; u 0 , β) of (1)-(2) depends continuously on (u 0 , β). The proof of this fact proceeds as follows: By exploiting the triangle inequality and the stability estimates (23)  and we are done. Here u 0 andũ 0 denote two initial conditions of (1), whileβ and β are two steepness parameters. The function g(T ) is increasing, and we conclude that the mapping where U is defined in (26)-(27), is continuous.