Regularization of ill-posed point neuron models

TL;DR

This paper investigates how regularizing steep firing rate functions in point neuron models can lead to well-posed equations and examines the convergence of these solutions to the ill-posed limit as the steepness increases, highlighting conditions for convergence.

Contribution

It provides a mathematical analysis of regularization in point neuron models, demonstrating convergence properties and conditions for the limit problem with a Heaviside firing rate.

Findings

Regularized models are well-posed and solvable with finite precision.

Solutions of regularized models converge subsequentially to the limit.

Convergence depends on measure-zero conditions for firing threshold crossing.

Abstract

Point neuron models with a Heaviside firing rate function can be ill-posed. That is, the initial-condition-to-solution map might become discontinuous in finite time. If a Lipschitz continuous, but steep, firing rate function is employed, then standard ODE theory implies that such models are well-posed and can thus, approximately, be solved with finite precision arithmetic. We investigate whether the solution of this well-posed model converges to a solution of the ill-posed limit problem as the steepness parameter, of the firing rate function, tends to infinity. Our argument employs the Arzel\`{a}-Ascoli theorem and also yields the existence of a solution of the limit problem. However, we only obtain convergence of a subsequence of the regularized solutions. This is consistent with the fact that we show that models with a Heaviside firing rate function can have several solutions. Our analysis assumes that the Lebesgue measure of the time the limit function, provided by the Arzel\`{a}-Ascoli theorem, equals the threshold value for firing, is zero. If this assumption does not hold, we argue that the regularized solutions may not converge to a solution of the limit problem with a Heaviside firing function.

Equations

$${\tau }_i{u}_i^{\prime }(t)=-u_i(t)+\sum _{j=1}^N{\omega }_{i,j}{S}_{\beta }[u_j(t)-u_{\theta }]+q_i(t),t\in (0,T],$$

$$u_i(0)=u_{\mathrm{init},i},$$

$$u_i(t)\in \mathrm{I}\mathrm{R},t\in [0,T],i=1,2,\dots ,N,$$

$$q_i(t)\in \mathrm{I}\mathrm{R},t\in (0,T],i=1,2,\dots ,N,$$

$$u_{\mathrm{init},i}\in \mathrm{I}\mathrm{R},i=1,2,\dots ,N,$$

$$u_{\theta }\in \mathrm{I}\mathrm{R},$$

$${\omega }_{i,j}\in \mathrm{I}\mathrm{R},i,j=1,2,\dots ,N,$$

$${\tau }_i\in {\mathrm{I}\mathrm{R}}_{+},i=1,2,\dots ,N,$$

$$\beta =1,2,\dots ,\infty ,$$

$$S_{\beta }[x]\text{mbox}isanapproximationoftheHeavisidefunctionH[x],$$

$$S_{\infty }[x]=H[x].$$

$$u_i(t)=u_{\beta ,i}(t),i=1,2,\dots ,N,$$

$${\mathbf{u}}_{\beta }(t)={\left(u_{\beta ,1}(t),u_{\beta ,2}(t),\dots ,u_{\beta ,N}(t)\right)}^T.$$

$$\mathbf{v}(t)=(v_1(t),v_2(t),\dots ,v_N(t){)}^T,$$

$$\tau {\mathbf{u}}^{\prime }(t)$$

$$\mathbf{u}(0)$$

$$\mathbf{u}(t)={\mathbf{u}}_{\beta }(t)\in {\mathrm{I}\mathrm{R}}^N,t\in [0,T],\text{mbox}see(\text{ref}u-beta),$$

$$\mathbf{q}(t)={\left(q_1(t),q_2(t),\dots ,q_N(t)\right)}^T\in {\mathrm{I}\mathrm{R}}^N,t\in (0,T],$$

$${\mathbf{u}}_{\theta }={\left(u_{\theta },u_{\theta },\dots ,u_{\theta }\right)}^T\in {\mathrm{I}\mathrm{R}}^N,$$

$${\mathbf{u}}_{\mathrm{init}}={\left(u_{\mathrm{init},1},u_{\mathrm{init},2},\dots ,u_{\mathrm{init},N}\right)}^T\in {\mathrm{I}\mathrm{R}}^N,$$

$$\omega =[{\omega }_{i,j}]\in {\mathrm{I}\mathrm{R}}^{N\times N},$$

$$\tau =\mathrm{diag}({\tau }_1,{\tau }_2,\dots ,{\tau }_N)\in {\mathrm{I}\mathrm{R}}^{N\times N}\text{mbox}isdiagonal,$$

$$\displaystyle S_{\beta}[\mathbf{x}]=(S_{\beta}[x_{1}],\ldots,S_{\beta}[x_{N}])^{T},\,\mathbf{x}={\color[rgb]{0,0,0}(x_{1},\ldots,x_{N})^{T}}\in\mathrm{I\!R}^{N}.$$

$$\Vert \mathbf{x}{\Vert }_{\infty }=\underset{1\le i\le N}{\max }|x_i|,\mathbf{x}=(x_1,...,x_N)\in {\mathrm{I}\mathrm{R}}^N.$$

$$|S_{\beta }(x)|<ϵ\text{mbox}forx<-\delta \text{mbox}and\beta >Q,$$

$$|1-S_{\beta }(x)|<ϵ\text{mbox}forx>\delta \text{mbox}and\beta >Q.$$

$$S(x)=\frac{1}{2}(1+\tanh (x)),$$

$$S_{\beta }[x]=S(\beta x).$$

$$\underset{\beta \in \mathrm{I}\mathrm{N},t\in [0,T]}{\sup }\Vert {\mathbf{q}}_{\beta }(t){\Vert }_{\infty }\le B<\infty ,B\in \mathrm{I}\mathrm{R},$$

$$\underset{\beta \to \infty }{\lim }{\mathbf{q}}_{\beta }(t)={\mathbf{q}}_{\infty }(t),t\in [0,T],$$

$$\underset{\beta \to \infty }{\lim }{\int }_0^t{\mathbf{q}}_{\beta }(s)ds={\int }_0^t{\mathbf{q}}_{\infty }(s)ds,t\in [0,T].$$

$$u_i^{\prime }(s)+{\tau }_i^{-1}{u}_i(s)={\tau }_i^{-1}\sum _{j=1}^N{\omega }_{i,j}{S}_{\beta }[u_j(s)-u_{\theta }]+{\tau }_i^{-1}{q}_i(s)$$

$${\left[u_i(s)e^{{\tau }_i^{-1}s}\right]}^{\prime }=e^{{\tau }_i^{-1}s}{\tau }_i^{-1}\sum _{j=1}^N{\omega }_{i,j}{S}_{\beta }[u_j(s)-u_{\theta }]+e^{{\tau }_i^{-1}s}{\tau }_i^{-1}{q}_i(s)$$

$$u_i(t)e^{{\tau }_i^{-1}t}=u_i(0)+{\int }_0^t{e}^{{\tau }_i^{-1}s}{\tau }_i^{-1}\sum _{j=1}^N{\omega }_{i,j}{S}_{\beta }[u_j(s)-u_{\theta }]ds+{\int }_0^t{e}^{{\tau }_i^{-1}s}{\tau }_i^{-1}{q}_i(s)ds.$$

$$|u_i(t)|e^{{\tau }_i^{-1}t}$$

$$=|u_i(0)|+\left(\sum _{j=1}^N|{\omega }_{i,j}|+\underset{s\in [0,T]}{\sup }|q_i(s)|\right)\left(e^{{\tau }_i^{-1}t}-1\right)$$

$$\le |u_i(0)|+\left(\sum _{j=1}^N|{\omega }_{i,j}|+B\right)\left(e^{{\tau }_i^{-1}t}-1\right),t\in (0,T],$$

$$\Vert {\mathbf{u}}_{\beta }(t){\Vert }_{\infty }\le \Vert {\mathbf{u}}_{\mathrm{init}}{\Vert }_{\infty }+\underset{i}{\max }\left(\sum _{j=1}^N|{\omega }_{i,j}|\right)+B,t\in [0,T].$$

$$\mathbf{v}=\underset{k\to \infty }{\lim }{\mathbf{u}}_{{\beta }_k}.$$

$$m(s;\mathbf{z})=\underset{j\in \{1,2,\dots ,N\}}{\min }|z_j(s)-u_{\theta }|,s\in [0,T].$$

$$Z(\mathbf{z})=\{s\in [0,T]|m(s;\mathbf{z})=0\}$$

Full text

Keywords: Point neuron models, ill-posed, regularization, existence.

1 Introduction

In this paper we analyze some mathematical properties of the following classical point neuron model:

[TABLE]

for $i=1,2,\ldots,N$, where

[TABLE]

Here, $u_{i}(t)$ represents the unknown electrical potential of the $i$th unit in a network of $N$ units. The nonlinear function $S_{\beta}$ is called the firing rate function, $\beta$ is the steepness parameter of $S_{\beta}$, $u_{\theta}$ is the threshold value for firing, $\{\omega_{ij}\}$ are the connectivities, $\{\tau_{i}\}$ are membrane time constants and $\{q_{i}(t)\}$ model the external drive/external sources, see, e.g., [2, 4, 5] for further details.

In computational neuroscience one often employs a steep sigmoid, or Heaviside, firing rate function $S_{\beta}$. This is due to both electrophysiological properties and mathematical convenience111Amari [1] analyzed the stationary solutions of neural field equations when $\beta=\infty$. . Unfortunately, the initial-condition-to-solution map for (1)-(2) can become discontinuous, in finite time, if a Heaviside firing rate function is used [7]. Such models are thus virtually impossible to solve with finite precision arithmetic [3, 15]. Also, in the steep, but Lipschitz continuous, firing rate regime, the error amplification can be extreme, even though a minor perturbation of the initial condition does not change which neurons that fire. It is important to note that this ill-posed nature of the model is a fundamentally different mathematical property than the possible existence of unstable equilibria, which typically also occur if a firing rate function with moderate steepness is used, see [7] for further details.

The solution of (1)-(2) depends on the steepness parameter $\beta$. That is,

[TABLE]

and the purpose of this paper is to analyze the limit process $\beta\rightarrow\infty$. This investigation is motivated by the fact that the stable numerical solution of an ill-posed problem is very difficult, if not to say impossible, see, e.g., [3, 15]. Consequently, such models must be regularized to obtain a sequence of well-posed equations which, at least in principle, can be approximately solved by a computer. Also, steep firing rate functions, or even the Heaviside function, are often used in simulations. It is thus necessary to explore whether the limit process $\beta\rightarrow\infty$ is mathematically sound. Similar type of studies, using different techniques, are presented in [8, 9] for the stationary solutions of neural field models.

In sections 3 and 4 we use the Arzelà-Ascoli theorem to analyze the properties of the sequence $\{\mathbf{u}_{\beta}\}$, where

[TABLE]

More specifically, we prove that this sequence has at least one subsequence which converges uniformly to a limit

[TABLE]

and that this limit satisfies the integral/Volterra version of (1)-(2) with $S_{\beta}=S_{\infty}$, provided that the Lebesgue measure of the time one, or more, of the component functions of $\mathbf{v}$ equals the threshold value $u_{\theta}$ for firing, is zero. Furthermore, in section 7 we argue that, if $\mathbf{v}$ does not satisfy this threshold property, then this function will not necessarily solve the limit problem.

According to the Picard-Lindelöf theorem [14], (1)-(2) has a unique solution, provided that $\beta<\infty$ and that the assumptions presented in the next section hold. In section 5 we show that this uniqueness feature is not necessary inherited by the limit problem obtained by employing a Heaviside firing rate function. It actually turns out that different subsequence of $\{\mathbf{u}_{\beta}\}$ can converge to different solutions of (1)-(2) with $S_{\beta}=S_{\infty}$. This is explained in section 6, which also contains a result addressing the convergence of the entire sequence $\{\mathbf{u}_{\beta}\}$.

For the sake of easy notation, we will sometimes write (1)-(2) in the form

[TABLE]

where

[TABLE]

2 Assumptions

Throughout this text we use the standard notation

[TABLE]

Concerning the sequence $\{S_{\beta}\}$ of finite steepness firing rate functions, we make the following assumption.

**Assumption **A

We assume that

a)

$S_{\beta}$, $\beta\in\mathrm{I\!N}$, is Lipschitz continuous,

b)

$0\leq S_{\beta}(x)\leq 1$, $x\in\mathrm{I\!R},\,\beta\in\mathrm{I\!N}$,

c)

for every pair of positive numbers $(\epsilon,\,\delta)$ there exists $Q\in\mathrm{I\!N}$ such that

[TABLE]

Reasonable/sound approximations of the Heaviside function satisfy A. For example, if $S_{\beta}$ is nondecreasing (for every $\beta\in\mathrm{I\!N}$), a) and b) hold and $\{S_{\beta}\}$ converges pointwise to the Heaviside function, then A holds. Also, if assumption A is satisfied and $\lim_{\beta\rightarrow\infty}S_{\beta}(0)=S_{\infty}(0)=H(0)$, then $\{S_{\beta}\}$ converges pointwise to the Heaviside function. Many continuous sigmoid approximations of the Heaviside function obey A. For example,

[TABLE]

We will consider a slightly more general version of the model than (4)-(5). More specifically, we allow the source term to depend on the steepness parameter, $\mathbf{q}=\mathbf{q_{\beta}}$, but in such a way that the following assumption holds.

**Assumption **B

We assume that $\mathbf{q}_{\beta}(t)$, $t\in[0,T]$, $\beta\in\mathrm{I\!N}\cup\{\infty\}$ is continuous and that

[TABLE]

Allowing the external drive to depend on the steepness parameter, makes it easier to construct illuminating examples. Please note that our theorems also will hold for the simplest case, i.e. when $\mathbf{q}$ does not change as $\beta$ increases.

In this paper we will assume that assumptions A and B are satisfied.

3 Uniformly bounded and equicontinuous

In order to apply the Arzelà-Ascoli theorem we must show that $\{\mathbf{u}_{\beta}\}$ constitute a family of uniformly bounded and equicontinuous functions. (For the sake of simple notation, we will write $u_{i}$ and $q_{i}$, instead of $u_{\beta,i}$ and $q_{\beta,i}$, for the component functions of $\mathbf{u}_{\beta}$ and $\mathbf{q}_{\beta}$, respectively). Multiplying

[TABLE]

with $e^{\tau_{i}^{-1}s}$ yields that

[TABLE]

and by integrating

[TABLE]

Hence, since $S_{\beta}[x]\in[0,1]$ and we assume that $\tau_{i}>0$ for $i=1,2,\ldots,N$,

[TABLE]

where the last inequality follows from (12). This implies that

[TABLE]

Since the right-hand-side of (15) is independent of $\beta$ and $t$, we conclude that the sequence $\{\mathbf{u}_{\beta}\}$ is uniformly bounded. Next, from the bound (15), the model equation (4), assumption (12) and the assumption that $S_{\beta}[x]\in[0,1]$ we find that also $\{\mathbf{u}^{\prime}_{\beta}\}$ is uniformly bounded. It therefore follows from the Mean Value theorem that $\{\mathbf{u}_{\beta}\}$ is a set of equicontinuous functions.

The Arzelà-Ascoli theorem [12, 6, 11] now asserts that there is a uniformly convergent subsequence $\{\mathbf{u}_{\beta_{k}}\}$:

[TABLE]

According to standard ODE theory, $\mathbf{u}_{\beta}$ is continuous for $\beta=1,2,\ldots<\infty$, and hence the uniform convergence implies that also $\mathbf{v}$ is continuous.

3.1 Threshold terminology

As we will see in subsequent sections, whether we can prove that $\mathbf{v}$ actually solves the limit problem with a Heaviside firing rate function, depends on $\mathbf{v}$’s threshold properties. The following concepts turn out to be useful.

For a vector-valued function $\mathbf{z}=(z_{1},z_{2},\ldots,z_{N})^{T}:[0,T]\rightarrow\mathrm{I\!R}^{N}$ we define

[TABLE]

Definition (Threshold simple)

A measurable vector-valued function $\mathbf{z}:[0,T]\rightarrow\mathrm{I\!R}^{N}$ is threshold simple if the Lebesgue measure of the set

[TABLE]

*is zero, i.e. $|Z(\mathbf{z})|=0$. *

Definition (Extra threshold simple)

A measurable vector-valued function $\mathbf{z}:[0,T]\rightarrow\mathrm{I\!R}^{N}$ is extra threshold simple if there exist open intervals

[TABLE]

such that

[TABLE]

With words, $\mathbf{z}$ is extra threshold simple if $|Z(\mathbf{z})|=0$ and the component functions of $\mathbf{z}$ only attains the threshold value for firing $u_{\theta}$ a finite number of times during $[0,T]$.

4 Convergence to the expected limit

4.1 Preparations

We will prove that the limit $\mathbf{v}$ in (16) solves the integral form of (4)-(5) with $S_{\infty}=H$, the Heaviside function, provided that $\mathbf{v}$ is threshold simple. The inhomogeneous nonlinear Volterra equation associated with (4)-(5) reads:

[TABLE]

where we consider the equations satisfied by the subsequence $\{\mathbf{u}_{\beta_{k}}\}$, see (16). We will analyze the convergence of the entire sequence in section 6. Note that we use the notation

[TABLE]

etc. in (19), see also (3) and (6).

The uniform convergence of $\{\mathbf{u}_{\beta_{k}}\}$ to $\mathbf{v}$ implies that the left-hand-side and the first term on the right-hand-side of (19) converge to the ”expected” limits as $k\rightarrow\infty$. Also, due to assumption (14), the third term on the right-hand-side does not require any extra attention. We will thus focus on the second term on the right-hand-side of (19).

Let, for $t\in[0,T]$ and $\delta>0$,

[TABLE]

where $m(s;\mathbf{v})$ is defined in (17), and $\mathbf{v}$ is the limit in (16). We note that, provided that $\delta>0$ is small, the set $r(\delta;t)$ contains the times where at least one of the components of $\mathbf{v}$ is close to the threshold value $u_{\theta}$ for firing. The following lemma turns out to be crucial for our analysis of the second term on the right-hand-side of (19)

Lemma 4.1

If the limit function $\mathbf{v}$ in (16) is threshold simple, then

[TABLE]

where $|r(\delta;t)|$ denotes the Lebesgue measure of the set $r(\delta;t)$.

Proof

• •

Since $\mathbf{v}$ is the uniform limit of a sequence of continuous functions, $\mathbf{v}$ is continuous and hence measurable.

• •

If $\mathbf{v}$ is threshold simple, then

[TABLE]

see (18).

• •

Let $t\in[0,T]$ be arbitrary.

• •

Assume that $\lim_{\delta\rightarrow 0^{+}}|r(\delta;t)|\neq 0$, or that this limit does not exist.

• •

Then $\exists\,\tilde{\epsilon}>0$ such that there is a sequence $\{\delta_{n}\}$ satisfying

[TABLE]

• •

By construction,

[TABLE]

and $|r(\delta_{1};t)|\leq T<\infty$. Hence,

[TABLE]

see, e.g., [11] (page 62). Since the sequence $\{|r(\delta_{n};t)|\}$ is nonincreasing and bounded below, $\lim_{n\rightarrow\infty}|r(\delta_{n};t)|$ exists.

• •

Next,

[TABLE]

i.e.

[TABLE]

• •

Hence,

[TABLE]

which contradicts (23).

4.2 Convergence of the integral

Lemma 4.2

If the limit $\mathbf{v}$ in (16) is threshold simple, then

[TABLE]

Proof

Let $t\in[0,T]$ and $\tilde{\epsilon}>0$ be arbitrary, and define

[TABLE]

From (22) we know that there exists $\Delta>0$ such that

[TABLE]

Choose a $\delta$ which satisfies $0<\delta<\Delta$. According to assumption A, for this $\delta$ and

[TABLE]

there exists $Q\in\mathrm{I\!N}$ such that (8) and (9) hold.

Recall that $\beta_{1},\,\beta_{2},\ldots,\,\beta_{k},\,\ldots$ are the values for the steepness parameter associated with the convergent subsequence $\{\mathbf{u}_{\beta_{k}}\}$ in (16). Let $K\in\mathrm{I\!N}$ be such that

[TABLE]

The existence of such a $K$ is assured by the uniform convergence of $\{\mathbf{u}_{\beta_{k}}\}$ to $\mathbf{v}$. From the definition of the set $p(2\delta;t)$, see (20) and (17),

[TABLE]

and from (28), and the triangle inequality, it follows that

[TABLE]

From (28) and (29) we find that

[TABLE]

Also, because of the properties of the Heaviside function,

[TABLE]

$j\in\{1,2,\ldots,N\}$. Consequently, due to (27) and assumption A, see (8) and (9), we find that

[TABLE]

Hence,

[TABLE]

for all $k>K$, where the second last inequality follows from (26), the fact that $|p(2\delta;t)|\leq T$ for $t\in[0,T]$ and (25). Since $\tilde{\epsilon}>0$ and $t\in[0,T]$ were arbitrary, we conclude that (24) must hold.

4.3 Limit problem

By employing the uniform convergence (16), the convergence of the integral (24) and assumption (14), we conclude from (19) that the limit function $\mathbf{v}$ satisfies

[TABLE]

provided that $\mathbf{v}$ is threshold simple. Recall that $\mathbf{v}$ is continuous. Consequently, if $\mathbf{v}$ is extra threshold simple, then it follows from the Fundamental Theorem of Calculus that $\mathbf{v}$ also satisfies the ODEs, except at time instances where one or more of the component functions equals the threshold value for firing:

[TABLE]

where $Z(\mathbf{v})$ is defined in (18).

The existence of a solution matter, for point neuron models with a Heaviside firing rate function, is summarized in the following theorem:

Theorem 4.3

If the limit $\mathbf{v}$ in (16) is threshold simple, then $\mathbf{v}$ solves (31). In the case that $\mathbf{v}$ is extra threshold simple, $\mathbf{v}$ also satisfies (32)-(33).

In [10] the existence issue for neural field equations with a Heaviside activation function is studied, but the analysis is different because a continuum model is considered. We would also like to mention that Theorem 4.3 can not be regarded as a simple consequence of Carathéodory’s existence theorem [13] because the right-hand-side of (32) is discontinuous with respect to $\mathbf{v}$.

5 Uniqueness

If $\beta<\infty$, then standard ODE theory [14] implies that (4)-(5) has a unique solution. Unfortunately, as will be demonstrated below, this desirable property is not necessarily inherited by the infinite steepness limit problem.

We will first explain why the uniqueness question is a subtle issue for point neuron models with a Heaviside firing rate function. Thereafter, additional requirements are introduced which ensure the uniqueness of an extra threshold simple solution.

5.1 Example: Several solutions

Let us study the problem

[TABLE]

where we assume that

[TABLE]

Note that the ODE (34) is not required to hold for $t=0$. Consider the functions

[TABLE]

Since

[TABLE]

it follows that both $v_{1}$ and $v_{2}$ solves (34)-(35).

Furthermore, with

[TABLE]

we actually obtain a third solution of (34)-(35). More specifically, the stationary solution

[TABLE]

We conclude that models with a Heaviside firing rate function can have several solutions – such problems can thus become ill-posed. (In [7] we showed that the initial-condition-to-solution map is not necessarily continuous for such problems, and that the error amplification ratio can become very large in the steep, but Lipschitz continuous, firing rate regime). Note that switching to the integral form (31) will not resolve the lack of uniqueness issue for the toy example considered in this subsection.

We also remark that:

• •

If we define $S_{\infty}(0)=1/2$, then neither $v_{1}$ nor $v_{2}$ satisfies the ODE (34) for $t=0$. (In the case $\omega=2u_{\theta}$, $v_{3}$ satisfies (34) for $t=0$.)

• •

If we define $S_{\infty}(0)=1$, then $v_{1}$, but not $v_{2}$, satisfies (34) also for $t=0$.

• •

If we define $S_{\infty}(0)=0$, then $v_{2}$, but not $v_{1}$, satisfies (34) also for $t=0$.

5.2 Enforcing uniqueness

In order to enforce uniqueness, we need to impose further restrictions. It turns out that it is sufficient to require that the derivative is continuous from the right and that the ODEs also must be satisfied whenever one, or more, of the component functions equals the threshold value for firing:

[TABLE]

Note that the ODEs (39) also must be satisfied for $t=0$, in case one of the components of $\mathbf{u}_{\mathrm{init}}$ equals $u_{\theta}$.

Definition (Right smooth)

A vector-valued function $\mathbf{z}:[0,T]\rightarrow\mathrm{I\!R}^{N}$ is right smooth if $\mathbf{z}^{\prime}$ is continuous from the right for all $t\in[0,T)$.

Theorem 5.1

Equations (39)-(40) can at the most have one solution which is both extra threshold simple and right smooth.

Proof

Let $\mathbf{v}$ and $\mathbf{\tilde{v}}$ be two solutions of (39)-(40) which are both right smooth and extra threshold simple:

[TABLE]

and

[TABLE]

where $I_{1},I_{2},\ldots,I_{L}$ and $\tilde{I}_{1},\tilde{I}_{2},\ldots,\tilde{I}_{\tilde{L}}$ are disjoint open intervals, see (17) and the definition of extra threshold simple in subsection 3.1.

Then there exist disjoint open intervals $\hat{I}_{1},\hat{I}_{2},\ldots,\hat{I}_{\hat{L}}$ such that

[TABLE]

Let us focus on one of these intervals, $\hat{I}_{l}=(a_{l},a_{l+1})$. Define

[TABLE]

and assume that

[TABLE]

which obviously holds for $l=1$. Then,

[TABLE]

where

[TABLE]

Note that, due to (41), $\mathbf{\gamma}(t)$ equals a constant vector $\mathbf{c}$, with components $-1,0$ or $1$, except possibly at $t=a_{l},\,a_{l+1}$:

[TABLE]

Furthermore, from (42) we find that

[TABLE]

Putting $t=a_{l}$ in (43) and invoking (44) and (46) yield that

[TABLE]

and from the right continuity of $\mathbf{d}^{\prime}$ and $\mathbf{d}$, (43), (44) and (45) we find that

[TABLE]

Since $\mathbf{\omega}\mathbf{\gamma}(t)=\mathbf{\omega}\mathbf{c}=\mathbf{0}$, $t\in(a_{l},a_{l+1})$, and $\mathbf{\omega}\mathbf{\gamma}(a_{l})=\mathbf{0}$, see (46), we conclude from (43)-(44) that $\mathbf{d}$ satisfies

[TABLE]

which has the unique solution $\mathbf{d}(t)=\mathbf{0}$, $t\in[a_{l},a_{l+1})$. Both $\mathbf{v}(t)$ and $\mathbf{\tilde{v}}(t)$ are differentiable on $[0,T]$ and hence continuous. It follows that, by employing the continuity of $\mathbf{v}$ and $\mathbf{\tilde{v}}$ at time $t=a_{l+1}$,

[TABLE]

Since $\mathbf{v}(a_{l+1})=\mathbf{\tilde{v}}(a_{l+1})$, we can repeat the argument on the next interval $[a_{l+1},a_{l+2}]$, and it follows by induction that $\mathbf{v}(t)=\mathbf{\tilde{v}}(t),\,t\in[0,T]$.

We would like to comment the findings presented in the bullet-points at the end of subsection 5.1 in view of Theorem 5.1: In order to enforce uniqueness for the solution of (34)-(35), we can require that the ODE (34) also should be satisfied for $t=0$. Nevertheless, this might force us to define $S_{\infty}(0)\neq\frac{1}{2}$, which differs from the standard definition of the Heaviside function $H$.

More generally, if one has accomplished to compute an extra threshold simple and right smooth function $\mathbf{v}$ which satisfies (31), then one can attempt to redefine $S_{\infty}[\mathbf{v}(t)-\mathbf{u_{\theta}}]$, $t\in\{a_{1},\,a_{2},\ldots,a_{L+1}\}$, such that (39)-(40) hold, and $\mathbf{v}$ is the only solution to this problem. This may imply that $S_{\infty}[\mathbf{v}(t)-\mathbf{u_{\theta}}]$ can not be generated by using the composition $H\circ[\mathbf{v}(t)-\mathbf{u_{\theta}}]$. Instead one must determine $z_{j,k}=S_{\infty}[v_{j}(a_{k})-u_{\theta}]$, $j=1,\,2,\ldots,\,N$, $k=1,\,2,\ldots,\,L+1$. More precisely, for each $k\in\{1,\,2,\ldots,\,L+1\}$ one gets a linear system of algebraic equations

[TABLE]

which will have a unique solution $(z_{1,k},\,z_{2,k},\,\ldots,\,z_{N,k})^{T}$ if the connectivity matrix $\mathbf{\omega}=[\omega_{i,j}]$ is nonsingular. (In this paragraph, $\{0=a_{1},\,a_{2},\ldots,a_{L+1}=T\}$ are the time instances one, or more, of the component functions of $\mathbf{v}$ potentially equals the threshold value for firing, see the definition of extra threshold simple in subsection 3.1).

6 Convergence of the entire sequence

We have seen that point neuron models with a Heaviside firing rate function can have several solutions. One therefore might wonder, can different subsequences of $\{\mathbf{u}_{\beta}\}$ converge to different solutions of the limit problem? In this section we present an example which shows that this can happen, even though the involved sigmoid functions satisfy assumption A.

6.1 Example: Different subsequences can converge to different solutions

Let us again consider the initial value problem (34)-(35), which we discussed in subsection 5.1. A finite steepness approximation of this problem, using the notation $u(t)=u_{\beta}(t)$, reads:

[TABLE]

where

[TABLE]

and $S_{\beta}$ is, e.g., either the $\tanh$-based sigmoid function (10)-(11) or

[TABLE]

Note that $\{\bar{S}_{\beta}\}$ converges pointwise, except for $x=0$, to the Heaviside function $H$ as $\beta\rightarrow\infty$. In fact, $\{\bar{S}_{\beta}\}$ satisfies assumption A.

We consider the case $\omega=2u_{\theta}$, and (34)-(35) therefore has three solutions $v_{1},\,v_{2}$ and $v_{3}$, see (36), (37) and (38) in subsection 5.1. Note that

[TABLE]

has the property

• •

$u^{\prime}_{\beta}(0)>c$ if $\beta$ is even,

• •

$u^{\prime}_{\beta}(0)<-c$ if $\beta$ is odd,

where $c>0$ is a constant which is independent of $\beta$. It therefore follows, argument not included, that

[TABLE]

and no subsequence converges to the third solution $v_{3}$. Figure 1 shows numerical solutions of (47)-(48) with steepness parameter $\beta=10\,000\,000,\\ \,10\,000\,001$, using the firing rate function (49) to define $\bar{S}_{\beta}$. (If one instead employs (10)-(11) in the implementation of $\bar{S}_{\beta}$, the plots, which are not included, are virtually unchanged).

We would like to mention that we have not been able to construct an example of this kind for Lipschitz continuous firing rate functions which converge pointwise to the Heaviside function also for $x=0$.

6.2 Entire sequence

We have seen that almost everywhere convergence of the sequence of firing rate functions to the Heaviside limit is not sufficient to guarantee that the entire sequence $\{u_{\beta}\}$ converges to the same solution of the limit problem. Nevertheless, one has the following result:

Theorem 6.1

Let $\mathbf{v}$ be the limit function in (16). If the limit of every convergent subsequence of $\{\mathbf{u}_{\beta}\}$ is extra threshold simple, right smooth and satisfies (39)-(40), then the entire sequence $\{\mathbf{u}_{\beta}\}$ converges uniformly to $\mathbf{v}$.

Proof

Suppose that the entire sequence $\{\mathbf{u}_{\beta}\}$ does not converge uniformly to $\mathbf{v}$. Then there is an $\epsilon>0$ such that, for every positive integer $M$, there must exist $\mathbf{u}_{\beta_{l}}$, $\beta_{l}>M$, satisfying

[TABLE]

The subsequence $\{\mathbf{u}_{\beta_{l}}\}$ can thus not converge uniformly to $\mathbf{v}$, but constitute a set of uniformly bounded and equicontinuous functions, see section 3. According to the Arzelà-Ascoli theorem, $\{\mathbf{u}_{\beta_{l}}\}$ therefore possesses a uniformly convergent subsequence $\{\mathbf{u}_{\beta_{l_{n}}}\}$,

[TABLE]

Due to (52),

[TABLE]

On the other hand, both $\mathbf{v}$ and $\mathbf{\tilde{v}}$ are limits of subsequences of $\{\mathbf{u}_{\beta}\}$, and are therefore by assumption extra threshold simple, right smooth and satisfies (39)-(40). Hence, Theorem 5.1 implies that $\mathbf{\tilde{v}}=\mathbf{v}$, which contradicts (53). We conclude that the entire sequence $\{\mathbf{u}_{\beta}\}$ must converge uniformly to $\mathbf{v}$.

One might argue that Theorem 6.1 only is of theoretical interest because it seems very difficult to guarantee that “the limit of every convergent subsequence of $\{\mathbf{u}_{\beta}\}$ is extra threshold simple, right smooth and satisfies (39)-(40)”.

7 Example: Threshold advanced limits

We will now show that threshold advanced limits, i.e. limits which are not threshold simple, may possess some peculiar properties. More precisely, such limits can potentially occur in (16), and they do not necessarily satisfy the limit problem obtained by using a Heaviside firing rate function.

With source terms which do not depend on the steepness parameter $\beta$, we have not managed to construct an example with a threshold advanced limit $\mathbf{v}$. If we allow $\mathbf{q}=\mathbf{q}_{\beta}$, this can, however, be accomplished as follows. Let

[TABLE]

where we, for the sake of simplicity, work with the firing rate function (49). Then,

[TABLE]

and we find that

[TABLE]

solves

[TABLE]

where

[TABLE]

It follows that

[TABLE]

and since, for any $\beta\in\mathrm{I\!N}$,

[TABLE]

we conclude that

[TABLE]

Note that

[TABLE]

but $\bar{v}(t)=u_{\theta}$ does not solve the limit problem

[TABLE]

because

[TABLE]

This argument assumes that $S_{\infty}[0]=1/2$. If one instead defines $S_{\infty}[0]=1$, then $\bar{v}$ would solve the limit problem.

Due to the properties of the firing rate function (49), the source term $q_{\beta}$ in (59) becomes discontinuous. This can be avoided by instead using the smooth version (10)-(11), but then the analysis of this example becomes much more involved.

The author does not know whether it is possible to impose restrictions which would guarantee that the limit $\mathbf{v}$ in (16) is threshold simple or extra threshold simple. This seems to be an herculean task.

8 Discussion and conclusions

If a Heaviside firing rate function is used, then the model (1)-(2) may not only have several solutions, but the initial-condition-to-solution map for this problem can become discontinuous [7]. It is thus virtually impossible to develop reliable numerical methods, which employ finite precision arithmetic, for such problems. One can try to overcome this issue by:

a)

Attempting to solve the ill-posed equation with symbolic computations.

b)

Regularize the problem.

As far as the author knows, present symbolic techniques are not able to handle strongly nonlinear equations of the kind (1), even when $\beta<\infty$. We therefore analyzed the approach b), using the straightforward regularization technique obtained by replacing the Heaviside firing rate function by a Lipschitz continuous mapping. This yields an equation which is within the scope of the Picard-Lindelöf theorem and standard stability estimates for ODEs. That is, well-posed and, at least in principle, approximately solvable by numerical methods.

Our results show that the sequence $\{\mathbf{u}_{\beta}\}$ of regularized solutions will have at least one convergent subsequence. The limit, $\mathbf{v}$, of this subsequence will satisfy the integral/Volterra form (31) of the limit problem, provided that the Lebesgue measure of the time one, or more, of the component functions of $\mathbf{v}$ equals the threshold value $u_{\theta}$ for firing, is zero. Unfortunately, it seems to be very difficult to impose restrictions which would guarantee that $\mathbf{v}$ obeys this threshold property, which we refer to as threshold simple. Also, the example presented in section 7 shows that, if the limit $\mathbf{v}$ is not threshold simple, then this function may not solve the associated equation with a Heaviside firing rate function.

One could propose to overcome the difficulties arising when $\beta=\infty$ by always working with finite slope firing rate functions. This would potentially yield a rather robust approach, provided that the entire sequence $\{\mathbf{u}_{\beta}\}$ converges, because increasing a large $\beta$ would still guarantee that $\mathbf{u}_{\beta}$ is close to the unique limit $\mathbf{v}$. However, the fact that different convergent subsequences of $\{\mathbf{u}_{\beta}\}$ can converge to different solutions of the limit problem, as discussed in section 6, suggests that this approach must be applied with great care. In addition, the error amplification in the steep firing rate regime can become extreme [7], and the accurate numerical solution of such models is thus challenging.

What are the practical consequences of our findings? As long as there does not exist very reliable biological information about the size of the steepness parameter $\beta$, and the shape of the firing rate function $S_{\beta}$, it seems that we have to be content with simulating with various $\beta<\infty$. If one observes that $\mathbf{u}_{\beta}$ approaches a threshold advanced limit, as $\beta$ increases, or that the entire sequence does not converge, the alarm bell should ring. All simulations with large $\beta$ must use error control methods which guarantee the accuracy of the numerical solution – we must keep in mind that we are trying to solve an almost ill-posed problem.

Competing interests

The author declares that he has no competing interests.

Acknowledgements

This work was supported by the The Research Council of Norway, project number 239070. The author would like to thank Prof. Wyller for several interesting discussions about the research presented in this paper.