Neural Field Models with Threshold Noise

The original neural field model of Wilson and Cowan is often interpreted as the averaged behaviour of a network of switch like neural elements with a distribution of switch thresholds, giving rise to the classic sigmoidal population firing-rate function so prevalent in large scale neuronal modelling. In this paper we explore the effects of such threshold noise without recourse to averaging and show that spatial correlations can have a strong effect on the behaviour of waves and patterns in continuum models. Moreover, for a prescribed spatial covariance function we explore the differences in behaviour that can emerge when the underlying stationary distribution is changed from Gaussian to non-Gaussian. For travelling front solutions, in a system with exponentially decaying spatial interactions, we make use of an interface approach to calculate the instantaneous wave speed analytically as a series expansion in the noise strength. From this we find that, for weak noise, the spatially averaged speed depends only on the choice of covariance function and not on the shape of the stationary distribution. For a system with a Mexican-hat spatial connectivity we further find that noise can induce localised bump solutions, and using an interface stability argument show that there can be multiple stable solution branches.


Introduction
The study of waves, bumps and patterns in models of Wilson-Cowan type [1] is now a very mature branch of mathematical neuroscience, as discussed in the review by Bressloff [2], with many practical applications to topics including working memory, visual processing, and attention. For a recent and comprehensive description of neural fields and their applications we refer the reader to the book [3]. It is only relatively recently that stochastic effects in neural fields have begun to be considered, with important applications to problems such as binocular rivalry waves [4] and perceptual switching [5]. These stochastic models are often obtained by considering the addition of noisy currents (notionally a "Gaussian random noise") to standard (deterministic) neural fields, and the resulting models are cast as stochastic nonlinear integro-differential equations driven by a Wiener process, such as in [6][7][8][9][10][11][12][13][14]. A rigorous probabilistic framework in which to study these equations has recently been provided by Faugeras and Inglis [15]. The analysis of patterns, waves and bumps in such models has been possible utilising tools from stochastic centre manifold theory (especially tools for weak noise analysis), Fokker-Planck reductions, and other techniques from stochastic calculus developed previously for PDEs. For a recent perspective on this approach the book by Bressloff is a highly valuable resource [16], as well as the paper by Inglis and MacLaurin [17], which presents a general framework in which to rigorously study the effect of spatio-temporal noise on travelling wave fronts. Indeed there is now a quite elegant body of rigorous theory growing up around neural field models with multiplicative stochastic forcing, as exemplified in the paper by Krüger and Stannat [18] using multiscale analysis, which moves beyond formal perturbation methods, to understand front propagation in particular. However, the original work of Wilson and Cowan suggests that another, perhaps more natural, way to introduce stochasticity into neural field models is by treating some of the system parameters as random variables. Indeed, threshold noise in a linear integrate-and-fire model is able to fit real firing patterns observed in the sensory periphery [19]. The simplicity of such models is also appealing from a theoretical perspective, and for a threshold described by an Ornstein-Uhlenbeck process it has recently been shown that analytical (and non-perturbative) expressions for the first-passage time distribution can be obtained [20].
To appreciate the original idea of Wilson and Cowan that threshold noise in switching networks can give rise to a probabilistic interpretation of network dynamics in terms of a smooth firing-rate function it is enough to consider a simple discrete time model for the evolution of neural activity x i (t), i = 1, . . . , N, in a network with connections w ij : Here H is a Heaviside step function with threshold h. If we now associate a threshold h i with each individual node and treat it as a random variable drawn from a normalised stationary probability distribution φ(h i ) at each time step then we can take the ensemble average of the above and find where Thus we obtain a smooth nonlinear deterministic model describing the average behaviour of a set of switch like elements with random thresholds, with the link between the two determined by the relationship f = φ. Since φ is a probability distribution, this relationship immediately implies a monotonically increasing firing-rate function. Given a realisation of the thresholds h i at some time t, it is of interest to ask how the spatial covariance structure of these random thresholds affects network dynamics. This is precisely the question we wish to address in this paper for continuum models of Wilson-Cowan type, in which the random firing threshold is now described as spatially continuous quenched disorder. Although we will restrict our attention to a Gaussian covariance function, we shall consider a broad class of stationary distributions, and present practical techniques from applied mathematics and statistics for working with non-Gaussian distributions. Moreover, by working with the Heaviside choice, as in (1), we will be able to build on the interface approach of Amari [21] to obtain explicit results for travelling fronts and bumps, and their dependence on the threshold noise structure.
In Sect. 2 we introduce our neural field model of choice, as well as the form of the stochastic threshold, namely its steady state distribution and spatial covariance structure. In Sect. 3 we show that, for a given realisation of the threshold, we may use the Amari interface approach to determine the instantaneous speed of a travelling front. We further show how to calculate the effects of the quenched spatial disorder arising from the noisy threshold using a perturbative approach, valid for small deviations of the threshold from its average value. We extend the approach for fronts to tackle stationary bumps in Sect. 4, where we also show how to determine the linear stability of localised solutions. This leads to a prediction that noise can induce multiple stable bumps, which we confirm numerically. Indeed throughout the paper we use direct numerical simulations to illustrate the accuracy of all theoretical predictions. Finally in Sect. 5 we discuss natural extensions of the work in this paper.

The Model
For mathematical convenience it is often easier to work with spatially continuous models rather than lattice models of the type described by (1). We consider a neural field u = u(x, t) ∈ R, x ∈ [0, L], t ∈ R + , whose dynamics is given by The kernel w represents the anatomical connectivity, and we have chosen to include the nonlinearity within the spatial integration, though activity based models with the nonlinearity outside the spatial integration may also be considered with the techniques described below (and are qualitatively similar in their behaviour). As it stands the model given by (3) is a standard Amari neural field model for the choice that h, the firing threshold, is a constant function. In this case the model is well known to support travelling waves, including fronts [22] and localised bump states in systems with a mixture of excitation and inhibition. For a review of such behaviour see [23], and for a recent overview of neural field modelling in general see [3].
In this paper we shall consider the case that h is a spatially random function. Given the wealth of mathematical knowledge for Gaussian disorder it would be highly convenient to make this choice for the threshold. However, this is a non-physiological convenience that we would prefer to avoid. Indeed it is very natural to expect threshold noise to be bounded and unlikely to be best described by a symmetric distribution. As such we will consider both Gaussian and non-Gaussian disorder and in particular skewed exponential distributions and distributions with compact support. We shall explicitly model the random firing threshold h(x) as where h 0 > 0 corresponds to the mean of the threshold, and g(x) denotes the quenched disorder with symmetric, bounded and non-negative spatial covariance function C(x, y). We shall fix this to be a Gaussian shape such that C(x, y) = C(|x − y|), with Here κ is the correlation length of the quenched disorder. Note that the variance of the threshold is given by 2 σ 2 . There exists a sequence of non-negative real numbers, λ m , m ≥ 1, which are eigenvalues of the covariance operator, associated with a sequence of eigenfunctions, e m , m ≥ 1, according to that form a complete orthonormal basis so that we may represent g(x) by its Karhunen-Loève decomposition [24][25][26] g(x) = ∞ m=1 λ m α m e m (x).
Here the α m are uncorrelated random variables with zero mean and unit variance, i.e. E(α m ) = 0 and E(α m α n ) = δ mn . The properties of the α m ensure that the Karhunen-Loève representation captures the first and second moment of g(x) exactly. The latter result follows from the fact that so that C(x, y) has the expected spectral representation. When the correlation length κ is much smaller than the domain size L, the Karhunen-Loève decomposition of g(x) for the Gaussian covariance function (5) with periodic boundary conditions can be very well approximated by [25] g(x) = ∞ m=0 β m λ m e (1) Here we have split the eigenfunctions e m (x) in (7) into two sets e 1 m (x) and e 2 m (x), which read , ω m = 2πm/L and e (1) for any n, m. To complete the model setup we need to specify the random coefficients β m and γ m . They are determined by the local distribution φ(g) of the quenched disorder g(x). If φ(g) is Gaussian, it suffices to choose the β m and γ m as uncorrelated univariate Gaussian random variables, namely E(β m ) = 0 = E(γ m ) and E(β m β n ) = δ mn = E(γ m γ n ). Indeed there is a large variety of methods to simulate Gaussian disorder including autoregressive-moving-averages [27], circular embedding [28] spectral representations [29] or the Karhunen-Loève decomposition [24][25][26]. However, if φ(g) is non-Gaussian, then β m and γ m are not described by a scaled version of φ(g). Thus we require suitable techniques to generate non-Gaussian disorder. Compared to Gaussian statistics, there are only a few methods for simulating non-Gaussian disorder. Amongst them, translation processes, in which a suitably chosen Gaussian model is non-linearly mapped to the desired non-Gaussian disorder, and a Karhunen-Loève decomposition feature most prominently [30,31]. We have chosen to employ a Karhunen-Loève decomposition as it is more robust [32], which means that we may use the same technique for both Gaussian and non-Gaussian disorder. We implement the method developed in [30] and provide details of the algorithm in Appendix A, though the main idea is as follows. As a starting point choose uncorrelated β m and γ m from the probability distribution φ(g) and generate a large number of samples of the quenched disorder g(x) according to (7). Since the β m and γ m are uncorrelated g(x) possesses the prescribed covariance function C(x, y). However, the probability distribution of g(x) differs from φ(g). This can be corrected by determining a new set of the β m and γ m , but these β m and γ m are now correlated. It is then possible to decorrelate the β m and γ m without changing their statistics, which in turn ensures that the statistics of g(x) still comply with φ(g). Overall, the method in [30] provides an iterative scheme such that the probability distribution of g(x) converges towards the prescribed distribution φ(g), while keeping the chosen covariance function exact in every iteration step. In Fig. 1 we show the three types of distribution that we use to realise threshold values. These are (i) a Gaussian distribution, (ii) a highly skewed shifted exponential distribution, and (iii) a piecewise linear distribution with compact support. The precise mathematical form for each of these is given in Appendix B.

Travelling Fronts
As mentioned in Sect. 1 much is now known about the effects of random forcing on neural field models. As regards travelling fronts the work of Bressloff and Webber [9] has shown that this can result in 'fast' perturbations of the front shape as well as a 'slow' horizontal displacement of the wave profile from its uniformly translating position. A separation of time-scales method is thus ideally suited to analysing this phenomenon, though we also note that more numerical techniques based upon stochastic freezing [33] could also be utilised. In this section we will explore the effects of quenched or 'frozen' threshold noise on the properties of a travelling wave, and in particular its speed. For a symmetric choice of synaptic kernel w(x) = w(|x|), which decays exponentially, the one-dimensional model (3) with a constant threshold is known to support a travelling front solution [22,23] that connects a high activity state to a low activity Equation (9)  state. In this case it is natural to define a pattern boundary as the interface between these two states. One way to distinguish between the high and the low activity state is by determining whether u is above or below the firing threshold. When denoting the position of the moving interface by x 0 (t), the above notion leads us to the defining equation Here, we are assuming that there is only one point on the interface as illustrated in Fig. 2, though in principle we could consider a set of points. For the choice (5) we see that C(x) is differentiable at x = 0, which means that the random threshold is differentiable in the mean square sense. The differentiation of (12) gives an exact expression for the velocity of the interface c in the form which modestly extends the original approach of Amari [21] with the inclusion of the term for h x . We can now describe the properties of a front solution solely in terms of the behaviour at the front edge that separates high activity from low, as described in [34,35]. To see this, let us consider a right moving front for which u( Hence, For simplicity we make the choice w(x) = exp(−|x|)/2 so that from (14) we find by differentiation with respect to x that for a right moving wave (for large t) By noting that and inserting (16) and (17) into (13), we find the wave speed c + of a right moving wave .
When we repeat the above derivation for a left moving wave, the wave speed c − is given by where we used Note that in the case of a constant threshold with which recovers a previous result, as discussed in [16]. If the front is moving to the right we have an exact expression for the speed (see also [36]): Examples of this relationship are shown in Figs. 3 and 4 where we plot both c(x) and the instantaneous front velocity extracted from a numerical simulation of (3). Figure 3 depicts results when the local probability distribution is a Gaussian for two different values of the correlation length κ, while Fig. 4 illustrates travelling fronts for thresholds that are locally distributed as a skewed exponential and a bump. We see excellent agreement between the numerical values and the expression (21).

Perturbative Calculation of Wave Speed
We can also perturbatively calculate the effects of threshold disorder on the speed of a travelling front. Substituting (4) into (21) and taking 1 we find that where the prime denotes differentiation. We will now take the spatial and ensemble average of (22). It is convenient to introduce an angle bracket notation to denote spatial averaging according to · ≡ L 0 ·dx/L. We find from (9) that g = β 0 √ λ 0 /L (since only the constant eigenfunction e (1) 0 (x) contributes to the integral and all other terms in (9) integrate to zero because of periodicity) and g = 0 (because of periodicity). Hence, the spatial average of c takes the compact form where L g 2 = ∞ m=0 β 2 m λ m + ∞ m=1 γ 2 m λ m and L g 2 = ∞ m=1 β 2 m λ m ω 2 m + ∞ m=1 γ 2 m λ m ω 2 m . Now taking expectations over the β m and γ m we obtain E( g ) = 0, LE( g 2 ) = λ 0 + 2 ∞ m=1 λ m and LE( g 2 ) = 2 ∞ m=1 λ m ω 2 m . Hence, we may construct c = E( c ) as This expression gives the lowest order correction term to the expression for speed (for a right moving wave) when the threshold is constant and takes the value h 0 . The term in square brackets in (24) is positive, and thus spatial disorder will always increase the average speed. This term also increases as the correlation length decreases since the λ m decay more slowly for smaller correlation lengths. Note, however, that the correction term remains uniformly small since λ m ∼ κ for all m when κ 1. Figure 5 shows c as a function of for a Gaussian distribution of threshold values, as well as results from directly averaging realisations of (21). Here the β m and γ m are randomly chosen from the unit normal distribution. We obtain identical results when these are uniformly distributed on [− √ 3, √ 3] (i.e. with mean zero and variance 1). The difference between using bounded distributions for the β m versus unbounded is that the maximum value of h is then bounded/unbounded. The results based on (24) are almost identical to those obtained from (21) for small values of , while minor deviations appear as we increase . In Fig. 6 we plot c as a function of for the shifted exponential and bump distribution. In addition, we show results for a Gaussian distri- bution with the same variance. We again observe very good agreement between the small noise expansion (24) and averaging (21) for small values of . In addition, the curves for the Gaussian threshold and for the non-Gaussian thresholds obtained from averaging (21) almost agree, while the expansion (24) yields identical results for both kinds of threshold noise. The latter is a direct consequence of the bi-orthogonality of the Karhunen-Loève expansion. Equation (24) only depends on the eigenvalues of the covariance function and not on the properties of the local distributions.

Stationary Bumps
Neural fields of Amari type are known to support spatially localised stationary bump patterns when the anatomical connectivity function has a Mexican-hat shape. In a one dimensional spatial model, and in the absence of noise, it is known that pairs of bumps exist for some sufficiently low value of a constant threshold and that only the wider of the two is stable [21]. For random forcing it is possible to observe noise-induced drifting activity of bump attractors, which can be described by an effective diffusion coefficient (using a small noise expansion) [11] or by an anomalous sub-diffusive process in the presence of long-range spatio-temporal correlations [37]. However, it is also known that spatial disorder can act to pin states to certain network locations by breaking the (continuous) translation symmetry of the system, as described in [38] for neural field models with heterogeneous anatomical connectivity patterns. It is the latter phenomenon that we are interested in here, especially since for a disordered threshold that breaks translational symmetry, it is not a priori obvious how many bump solutions exist and what their stability properties are.
A one-bump solution q(x) is characterised by a width , such that for two values x 1 and x 2 with = x 2 − x 1 we have q(x) ≥ h(x) for x 1 ≤ x ≤ x 2 . Using (3) a onebump solution therefore satisfies We can determine the linear stability of bumps by studying the (linearised) evolution of a perturbation v(x)e λt around the stationary bump q(x). We hence find from (3) where Q(x) = q(x) − h(x). Demanding that the perturbations at x 1,2 be non-trivial yields the spectral equation det(A − (1 + λ)I 2 ) = 0, where I 2 is the identity matrix in R 2×2 and The eigenvalues are then given by λ = λ ± : Note further that for h (x) = 0 we have |Q (x 1 )| = |Q (x 2 )| = |w(0) − w( )| and therefore λ − = 0 as expected from translation invariance.

Simple Heterogeneity
We first consider a simple form of heterogeneity to present the ideas, and then move to more general heterogeneity. Suppose w(x) = e −α(1−cos x) − Be −β(1−cos x) and h(x) = h 0 + cos x, and the domain is [0, 2π]. Then we have bumps which have their maximum at either 0 or π . Suppose the maximum is at 0 and x 1 = −a for 0 < a < π. Since we need h(x 1 ) = h(x 2 ) and the threshold is symmetric around 0, we immediately arrive at Choosing α = 5, B = 0.76, β = 3, h 0 = 0.05 and = 0 we have two bump widths, as shown in Fig. 7 (a = /2). We find that the larger root is stable and the other is unstable, and both have a zero eigenvalue as expected. Increasing from zero breaks the translational invariance of the system and we obtain 4 bumps for = 0.01 as shown in Fig. 8: the two that exist for the homogeneous case ( = 0), now centred around x = π , and similarly two centred around x = 0. (We no longer restrict to

General Heterogeneity
We keep w(x) = e −α(1−cos x) − Be −β(1−cos x) and now consider a general h(x) without symmetries. The problem of the existence of a bump is this: given h(x) and a value of x 1 , find a value of x 2 such that h( is given by (26). This will not happen generically but only at isolated points in the (x 1 , x 2 ) plane. Thus we need to solve the simultaneous equations and we only search for solutions which satisfy 0 < x 1 < 2π and x 1 < x 2 < x 1 + 2π.
We now apply this general concept to the stochastic threshold (9). For each realisation and set of parameter values we use Newton's method with 1000 randomly chosen initial values which satisfy (34). Out of these 1000 initial values, the number of distinct solutions of (32)-(33) that satisfy (34) is recorded, and their stability is determined as described above. We then check these solutions to verify that q(x) > h(x) only for Any that do not satisfy this inequality are discarded. Typical results for a Gaussian threshold are depicted in Fig. 9 for κ = 1 and = 0.01. We see that width values ( ≡ x 2 − x 1 ) are clustered around two values (those corresponding to the homogeneous system) and that the only stable ones are those with large widths (inherited from the homogeneous case). Two example bumps, one stable and the other unstable, from Fig. 9 are shown in Fig. 10. Space-time simulations of (3) using these two bumps as initial conditions are shown in Fig. 11. The stable bump is stationary, as expected, while the solution starting close to (but not exactly at) the unstable bump rapidly increases in size and then slowly moves towards a stable solution. In Fig. 12 we show the number of solutions of (32)-(33) as well as the number and fraction of stable solutions as a function of . While the number of solutions exhibits an increasing trend, the number of stable solutions decreases, leading to an overall decrease in the fraction of stable solutions. When lowering the correlation length five times (Fig. 13), we observe a similar behaviour. Note, however, that the number of solutions has increased significantly. When we fix and vary κ the number of solutions decays quickly as shown in Fig. 14. At the same time, the number of stable solutions remains almost constant, resulting in a strong increase of the fraction of stable solutions. Overall, we find that varying the amplitude of the threshold heterogeneity by changing more strongly affects solution numbers when κ is small, and that the value of κ has a significant effect on how many fixed points exist (and a lesser effect on the number of those which are stable).
A more detailed view on bump solutions with a Gaussian threshold is presented in Fig. 15, where we show the distribution of values as a function of and the proportion which are stable. For = 0 the probability distribution consists of two delta functions at the stable and unstable bump, respectively. As increases two branches of solutions emerge from the solutions in the homogeneous case, which widen for larger values of . Note that as in the homogeneous case, only large-width bumps are stable. When = 0 bumps only exist below a critical value of h 0 , and a branch of stable bumps coalesces with a branch of unstable bumps at this critical value when h 0 is varied. Figure 16 shows results for a Gaussian threshold when we change h 0 for = 0.002. We again observe two solution branches that only exist below a critical value of h 0 . Each solution branch is smeared out compared to the homogeneous limit, indicating a probability distribution that has a finite and non-zero width.
Using the shifted exponential distribution for the threshold we obtain the results plotted in Fig. 17. We again observe two solution branches that emerge from the solutions in the homogeneous case as increases, with the stable solutions confined to the upper branch. In contrast to the Gaussian case in Fig. 15 the two solution branches do not grow symmetrically around the values of the homogeneous case. This is a manifestation of the highly skewed character of the exponential distribution compared to the symmetric Gaussian distribution. The probability distribution of the widths also exhibits much more structure compared to the Gaussian case.
We know that for a Heaviside firing rate, Mexican-hat connectivity, and a constant firing threshold, multibump solutions cannot exist [39]. However, breaking the translational invariance by using a heterogeneous threshold does allow such solutions to exist. We show an example of this behaviour in Fig. 18, where the L 2 norms of stable steady states of (3) are shown as a function of h 0 . Each point corresponds to a different realisation of h(x) and a different initial condition. We clearly see different "bands", each identified with an integer number of bumps in a solution, and the L 2 norm of these solutions increases as the number of bumps increases, in the same way as seen for systems with a smooth firing-rate function and homogeneous threshold [40][41][42], or with an oscillatory coupling function [39].

Conclusion
In this paper we have explored the role of threshold noise on travelling fronts and bumps in a simple neural field model with a Heaviside nonlinearity. For a frozen form of disorder and a given realisation of a spatial threshold the standard rules of calculus apply and we have exploited the interface approach of Amari to obtain exact results about solution properties for both existence and stability. The theory that we have developed is not restricted to any special choice of distribution for describing the threshold noise, apart from that sample trajectories be differentiable in the mean square sense. It is worth noting that the stochastic threshold model presented here is formally equivalent to already published stochastic nonlinear integrodifferential equations with constant threshold [6][7][8][9][10][11][12][13][14] when we employ the transfor- mation v(x, t) = u(x, t) − g(x). However, our approach permits the analysis of strong noise (see e.g. [20]) and hence will allow us to move beyond perturbative expansions.
Theoretical predictions have been shown to be in excellent agreement with numerical simulations of both Gaussian and non-Gaussian threshold models. As such we have a viable mathematically tractable model of a noisy neural tissue that would be of interest to explore in a variety of more neurobiologically relevant scenarios. For example, temporal correlations in excitability of neural tissue would be expected to strongly affect wave propagation and could be easily modelled with an appropriate choice for h = h(x, t). To investigate the consequences for wave speed all that would be required would be a minimal extension of the formula for interface dynamics (13) with the replacement of the numerator u t according to u t → u t − h t . Moreover, the inclusion of a linear adaptation current, as is common for describing spike frequency adaptation, would allow the study of travelling pulses as well as fronts, building on the deterministic approach in [43] and more recently in [44] for understanding cortical waves in epilepsy. The work here can also be extended to planar systems, allowing the study of spiral waves [45] and an investigation of how noise levels could be used to control the scale and size of a spiral, as suggested in [46]. These are topics of current investigation and will be reported upon elsewhere.
In the following it is convenient to introduce the notation α (k,l) m , which denotes the mth expansion coefficient at the kth iteration for the lth realisation of the quenched disorder. To initialise the scheme, we generate M sets of uncorrelated α (0,l) m , l = 1, . . . , M, drawn from the prescribed non-Gaussian distribution φ(g) shifted to mean zero and scaled to unit variance. A convenient way for doing this is to use inverse transform sampling based on the cumulative distribution function F φ of φ(g). Strictly speaking the α (0,l) m could be generated from any mean zero and unit variance distribution, but starting with the prescribed distribution φ(g) might speed up convergence. For numerical purposes, the sum in (7) will only contain N terms. One method to determine N is to impose the condition  i.e. that the maximal L 1 norm of the difference between the given covariance function and its approximation is smaller than some ε > 0. Once N is fixed, each iteration step proceeds as follows: 1. Generate M samples of the quenched disorder where g (k,l) denotes the lth realisation of the quenched disorder at the kth iteration. Note that because the α (k,l) m are uncorrelated with unit variance, the numerical error in determining the covariance function only depends on N (to satisfy (8)) and on M (for high-fidelity averaging). 2. Numerically determine the cumulative distribution function The mean of α (k+1,l) m is zero by construction, but the variance is unequal from one due to (39). Is therefore necessary to scale the α (k+1,l) m to have unit variance. 5. While the α (k+1,l) m now give rise to the prescribed probability distribution φ(g), they are correlated due to the mapping in (39). To decorrelate them, we use an iterative product-moment orthogonalisation technique.

Appendix B: Local Probability Distributions
We here provide details of the three zero-mean local probability distributions used in this study. We consider a Gaussian distribution a highly skewed shifted exponential distribution We refer to the last distribution as bump distribution. By demanding that it is normalised, we find α = 1/(L 2 − B 2 ). The variances of the exponential and bump distribution are given, respectively, by