Kernel Reconstruction for Delayed Neural Field Equations

Understanding the neural field activity for realistic living systems is a challenging task in contemporary neuroscience. Neural fields have been studied and developed theoretically and numerically with considerable success over the past four decades. However, to make effective use of such models, we need to identify their constituents in practical systems. This includes the determination of model parameters and in particular the reconstruction of the underlying effective connectivity in biological tissues. In this work, we provide an integral equation approach to the reconstruction of the neural connectivity in the case where the neural activity is governed by a delay neural field equation. As preparation, we study the solution of the direct problem based on the Banach fixed-point theorem. Then we reformulate the inverse problem into a family of integral equations of the first kind. This equation will be vector valued when several neural activity trajectories are taken as input for the inverse problem. We employ spectral regularization techniques for its stable solution. A sensitivity analysis of the regularized kernel reconstruction with respect to the input signal u is carried out, investigating the Fréchet differentiability of the kernel with respect to the signal. Finally, we use numerical examples to show the feasibility of the approach for kernel reconstruction, including numerical sensitivity tests, which show that the integral equation approach is a very stable and promising approach for practical computational neuroscience.


Introduction
In recent years, studying the activity of neural tissue and development of mathematical and numerical techniques to understand neural processes has led to improved neural field models. Since the early work of Wilson, Cowan and Amari in the 1970s neural field models have become an effective tool in neuroscience [1][2][3][4].
The neural networks occurring in nature are typically complex systems sporting a large variety of properties in space and time. Simplifying their analysis is generally difficult-in particular when one considers the many billions of neurons of the entire human nervous system, where each of these neurons can be considered as a complex biological system by and in itself, cf. [5]. However, neural field models describe these complicated system mathematically in a few equations, essentially by using the large number of neurons to achieve simplification in terms of mass action. Thus these models consider averages of the neural activity as a dynamical variable, and averages of neural properties as parameters. The derivation of neural models from properties of single neurons and their networks, and the analysis of the resulting activity, remains a major focus of current research [1,[6][7][8][9][10][11].
In this century, there are many papers on the neural field equation with and without delays. Some of the studies provide a framework for the existence, uniqueness and stability of the solutions of the neural field equation such as [8][9][10][11][12][13][14][15], while others consider building effective methods to investigate and assimilate the neural field activities, see for example [16][17][18][19][20][21] with techniques of Data assimilation and Inverse Problems applied to the case without delays. Recently, Nogaret et al. [7] built a model construction method using an optimization technique to assimilate neural data to determine parameters in a detailed neural model including delay.
A challenge often encountered in the study of living systems is to estimate a spatial connectivity kernel w. In a neural system this connectivity kernel usually corresponds to the synaptic footprint, i.e., the connections from a neuron to others by synapses forming between its branching axon and their dendritic trees. Typically, measurements are available for the activity function u at particular spatial locations, e.g., where neurons are patch-clamped or electrodes are placed in the extracellular medium. The task then becomes to derive the spatial connectivity from these experimental data. This approach limits the estimation of connectivity to the set of spatial locations of measurements. In the present work, we propose to improve this conventional approach by studying the inverse problem where the full activity function u is given at each location in a given spatial domain and the underlying spatial connectivity is derived. The problem of having limited measurements is part of subsequent work combining inverse techniques with state estimation techniques. Here, we focus on the problem to reconstruct the kernel w when u is known.
The present work considers neural field models that involve delayed spatial interactions and where the delay may depend on the distance between spatial locations [11,14,22]. We will assume that the delay function D(r, r ) between spatial locations r, r is known. For instance, this is the case when the delay is linked to the geometry of the problem, e.g., when D(r, r ) ∼ r − r , the distance between the points r and r in some domain Ω. This assumption is common in practice, since for direct neural connections the delay is essentially the distance divided by the signal propagation speed, which can be assumed to be a universal constant in a first approximation.
Neural field models consider spatially nonlocal interactions, which may be expressed equivalently either by higher orders of spatial derivatives or by spatial integrals [22,23]. In the first part of this paper, we will show how the methods used in [12] can be modified to study existence and stability of solutions in a neural field model with delay. The basic idea is to split the integral operators under consideration into parts with positive and negative temporal arguments. As a result we obtain a direct and flexible basic existence proof for a delay neural field equation, which includes a constructive method based on integral equations only. These results have been derived by other authors [8,10,11,24] with more sophisticated techniques, but it is non-trivial that the arguments used for neural fields without delay are applicable to the delay case, and the approach in our Sect. 2, based on several relatively simple functional analytic arguments, is of interest by itself.
Second, we will show that the kernel reconstruction problem for the delay neural field equation can be reformulated into a family of integral equations of the first kind. When several trajectories of neural activity are given, the family of integral equations is vector valued. This turns out to be an ill-posed problem, for smooth neural activity it is even exponentially ill posed. To formulate stable numerical methods for its solution, we need to employ regularization. Here, we use a spectral approach to classical Tikhonov regularization [25][26][27]. We then study the sensitivity of the mapping u → w showing that its regularized version is Fréchet differentiable, and we calculate the derivative by means of integral equations.
In the third part of the paper, we show by a numerical study that the kernel reconstruction from a delay neural field is feasible. We numerically solve the family of integral equations under consideration by a collocation method and provide a study of reconstructions based on the regularization of the ill-posed integral operators under consideration. This includes a study of the influence of measurement noise on the reconstruction quality and tests of the role of the regularization parameter.
We start with a concise version of the equations in Sect. 2, and in Sect. 3 prepare our inverse approach by a study of the existence for the delay neural field equation. The central section, Sect. 4, serves to develop a family of integral equations to solve the inverse problem for the delay neural field equation. The numerical realization of the approach is shown in Sect. 5, where we demonstrate that with an appropriate regularization the inverse problem is solvable, i.e., prescribed kernels can be constructed and reconstructed kernels generate a neural environment leading to the prescribed neural behaviour.

The Mathematical Model
In neural dynamics, neurons send electrical spikes to each other through axons terminating in synapses. Let u(r j , t) denotes the average membrane potential of the j th neuron located at position r j at time t in a network of N units. Let W (r j , r i ) be the average connectivity strength between neuron at position r i and neuron at position r j . The function f is the activation rate or firing rate function, which describes the conversion of the membrane potential u(r j , t) into a spike train S(r i , t) = f [u(r i , t)], which is then leading to an excitation of neurons at location r j with strength W (r j , r i )S(r i , t). The dynamics of the excitation is now described by the ODEs This combination of an exponential decay with characteristic time τ and a sum of excitation terms is commonly called a 'leaky integrator model'. The sum represents the net-input to unit j , i.e., the weighted sum of activity delivered by all units i that are connected to unit j with a connection strength W (r j , r i ); cf. [12,28]. The continuous version of (1) is obtained by considering neurons which are continuously distributed over the space Ω, e.g., in a plane with Ω ∈ R 2 or Ω ∈ R 3 and by replacing the sum by an integral. This leads to the simplest form of the Amari neural field equation [4], Here u(r, t) indicates a neural field representing the activity of the population of neurons at position r and time t. The second term on the right-hand side represents the synaptic input, where f is the activation (or firing rate) function of a single neuron. The kernel w(r, r ) is often referred to as the synaptic footprint [29][30][31] or the connectivity function [12,14,32,33]. It presents the strength of the connection between neurons located at r and r . The function w incorporates three different kinds of meaning: the existence of a connection in the first place, if w = 0, the functional effect of either excitation, if w > 0, or inhibition, if w < 0, and finally the strength of the connectivity via |w| [4,12,34]. Although the neural field equation (2) represents several biological mechanisms, this form still neglects any delay between spatial locations. In reality, finite transmission speeds in axons, synapses and dendrites cause a functionally significant delay. Taking it into account, the neural field equation involving delayed interactions becomes where the delay is typically assumed to be D(r, r ) D (r, r )/v, i.e., the total length of the neural fibersD connecting locations r and r , divided by v, the finite transmission speed of neural signals (action/post-synaptic potentials) along those fibers. In general, D is not constant but continuous. Equation (3) is accompanied by initial conditions. These depend on the geometry of the spatial domain and the specific temporal dynamics under study. They are considered in detail in the subsequent sections.
The existence of solutions to the neural field equation (3) has been investigated in various papers already [10][11][12]. For example, Potthast and beim Graben [12] provide the proof of existence and its analysis in the case of no delay, i.e. for D(r, r ) = 0. In addition, Faugeras and Faye [10], in their Theorem 3.2.1, state the general existence of solutions with a reference to the generic theory of delay equations, based on work such as [35]. We also point out the work of Van Gils et al. [8] employing the sun-star calculus for their analysis and [24] in which the local bifurcation theory for delayed neural fields was developed. Here, we develop arguments on how to use the basic functional analytic calculus to work for the delay case as well, with the goal to present a short and elementary approach which is easily accessible.

The Delay Neural Field Equation
In this work, we study the neural field equation (3) on some bounded domain Ω ⊂ R m in a space with dimension m = 2 or m = 3. We assume that the transmission delay D(r, r ) of neural excitation or inhibition between r and r is bounded on Ω × Ω, i.e. there is a constant c T such that At time t ∈ R, the neural fields u(r, t) at a point r ∈ Ω might receive excitations from the past with a maximal delay of c T . Working on the time interval [0, ρ] with ρ > 0, equation (3) is complemented by initial conditions in the time interval [−c T , 0]. The initial condition for the delay neural field equation is given by We lay ground for our inverse and sensitivity analysis by a basic derivation of the unique solvability of equation (3), using tools from functional analysis and integral equations. Our investigation here makes a smoothness assumption for the activity function f and the connectivity kernel w. We consider a continuous activation function f (s) for s ∈ R and an activation threshold η. This function may be interpreted as the mass action probability of neurons firing if their membrane potential is over the threshold, and can be derived from a stochastic neuron models [6,36]. Typically [1,29], f is approximated by the logistic sigmoidal function with some steepness parameter σ > 0 and threshold η. For the function f : R → R + we note that Here, we will work with general Lipschitz continuous functions f satisfying this condition. We assume that the kernel w satisfies such that we obtain a well-defined integral of the form The condition with some constant C 1 leads to g being bounded on Ω × R. We need g(r, s) to be continuous in dependence of r and s, which for continuous functions u and D is achieved by the additional condition Now, existence is given by the following result.
Proof We first need some preparations. We will need to split the function u(r, s − D(r, r )) into the part where the time variable . This is carried out by defining and χ − (r, t) := 1 − χ + (r, t). The function χ − is equal to 1 for negative time arguments and we have 1 = χ + + χ − . For studying the existence of solutions of the delay neural field equation (3) we define the operators and for r ∈ Ω and t ∈ [0, T ]. By integration with respect to time the solution of (3) can be reformulated as for r ∈ Ω and t ∈ [0, ρ] with an auxiliary parameter ρ. Differentiating equation (11) with respect to time, we return to the delay neural field equation (3). We can now split the operators as follows: where the last equality is obtained from Here, the function u(r, t) needs to be considered on Ω × [0, ρ] only and we can study the fixed-point equation in BC(Ω × [0, ρ]). Any solution to equation (13) will be continuously differentiable with respect to time and satisfy the delay neural field equation (3). We now show that for sufficiently small parameter ρ > 0 the operator A is a contraction on the space We will carry out these arguments in four steps, I-IV.
III. Integration of equation (18) with respect to t ∈ [0, ρ] leads to where · ρ as defined in equation (14). Now, for the operator A we obtain the estimate with In the case where ρ is small enough to guarantee that q < 1 by equation (20), we have shown that A is a contraction on BC(Ω × [0, ρ], · ρ ).
IV. According to the Banach fixed-point theorem, there is one and only one fixed point u * for the fixed-point equation (13). We have shown the existence of a unique solution u(x, t) for all t ∈ [0, ρ]. Now, the same argument applied to the interval [ρ, 2ρ] and subsequent intervals [2ρ, 3ρ] etc. in the same way. This leads to the existence and uniqueness result on the interval [0, T ].
Remark We note that the proof also works when some bounded continuous forcing term I (r, t), r ∈ Ω, t ∈ [0, T ], is added to the neural field equation (3). It leads to an additional term in Eq. (13), for which all arguments remain valid.
It is well known [21,27] that Banach's theorem also provides a constructive method to calculate the fixed point by successive iterations. Let u 1 be a starting function. Then the sequence defined by converges to the unique fixed point u * . An error estimate for this iteration process based on equation (20) is obtained from Induction immediately leads to the full error estimate For our numerical calculations we have, however, instead used Runge-Kutta or Euler methods applied to the differential form of the delay neural field equation.

The Inverse Problem of Kernel Reconstruction with Delays
We now come to the kernel reconstruction from given dynamical neural patterns with delay. We first formulate a regularized kernel reconstruction approach based on integral equations in Sect. 4.1, then we carry out a sensitivity analysis in Sect. 4.2.

Kernel Reconstruction with Delays
Usually, we will observe the dynamical evolution of some pattern for a system under consideration. More generally, observations may start from different inital patterns that lead to different dynamical trajectories in the phase space. If we have N such trajectories, the task is to find the kernel which will predict these trajectories when the N initial conditions are provided. In more detail, the goal of this section is to investigate the inverse problem of kernel reconstruction for the delay neural field equation (3). We assume that • the nonlinear activation function f : R → R + is known, and • the delay function D : Ω × Ω → [0, c T ] is given.
Here, we reformulate the inverse problem into a family of integral equations of the first kind and study their solution by regularization methods. As a first step, we define and for ξ = 1, 2, . . . , N. With the integral operator W defined by the inverse problem is reformulated as the equation with ξ = 1, 2, . . . , N, where the kernel w(r, r ) with r, r ∈ Ω of the linear operator W is unknown. Equation (28) can be written as with φ = (φ (1) , φ (2) , . . . , φ (N ) ) T and ψ = (ψ (1) , ψ (2) , . . . , ψ (N ) ) T , where we search for the unknown operator W . An alternative is to rewrite equation (28) as for every fixed r ∈ Ω with and w r := w r, r , r,r ∈ Ω.
Equation (30) is a family of integral equations for the unknown kernel w(r, r ), where each function w r = w(r, ·) provides a different integral equation with a different integral kernel and a different left-hand side. Its structure is given by the integral operator with kernel for r ∈ Ω. For N > 1 this kernel is a vector of functions φ (ξ ) (r , t − D(r, r )) with ξ = 1, . . . , N. Now, our inverse problem equation (30) is given by for r ∈ Ω. For each r ∈ Ω equation (35) is a Fredholm integral equation of the first kind with continuous kernel φ. The operator V r is a compact operator on the spaces . It is well known (cf. [21,25,27,37]) that this equation is ill posed, i.e. it does not need to have unique solutions and if it has a solution in general this solution does not depend continuously on the right-hand side. Ill-posed equations need some regularization method (cf. [26]) in order to obtain a stable solution. A standard approach to regularization is built on the singular system (cf. [27]) of the operator under consideration. In summary, for a compact linear operator A : X → Y between Hilbert spaces X and Y , and its adjoint A * , the singular values μ n of the operator A are the non-negative square roots of the eigenvalues of the self-adjoint compact operator A * A : X → X. This leads to a representation of the operator as a multiplication of two orthonormal systems g n : n ∈ N in X and y n : n ∈ N in Y . Hence, this corresponds to a spectral representation of the operator A in the form Ag = ∞ n=1 μ n g, g n y n , for g ∈ X. For the orthonormal systems g n and y n we obtain Ag n = μ n y n , A * y n = μ n g n .
Here, in the case A that is injective, the inverse of A is given by 1 μ n y, y n g n (38) or, if A is not injective, the inverse A −1 in equation (38) projects onto the orthogonal space N(A) ⊥ = {g| g, g * = 0, ∀g * ∈ N(A)}. Because of the compactness of the operators A, the singular values are a sequence mostly accumulating at zero. So, the behaviour of | 1 μ n | → ∞, n → ∞ enlarges small errors causing the instability of applying the inverse. The practical behaviour of the sequence of singular values μ n provides important insight into the nature of the instability. For the application at hand the problem is strongly ill posed for strong smoothness of the function φ.
To deal with this instability, we apply regularization techniques to minimize the value of the factor 1 μ n for large n. We replace it by another factor q n which is bounded for n ∈ N, and we modify the inverse operator by R α y = ∞ n=1 q (α) n y, y n g n , where α > 0 is known as regularization parameter and the specific choice of damping factors q (α) n := μ n α + μ 2 n , n∈ N (40) leads to the famous Tikhonov regularization (see for example [21,[25][26][27]). u(r, t) for r ∈ Ω and t ∈ [0, T ] be some neural activity function, which obeys the neural field equation (3) with true kernel w * and some initial conditions u(r, t) = u 0 (r, t) for (r, t) ∈ Ω × [−c T , 0]. Then the application of the Tikhonov regularization (39) to the integral equation (35) leads to the reconstruction w α (r, r ) of P w * , where P is applied to the second argument of w(r, r ) as the projection of w * r onto N(V r ) ⊥ , i.e., it is defined as

Theorem 4.1 Let
Proof Here, we base our reconstruction on a well-known result (cf. [21], Theorem 3. we see by w r = P w r + (I − P )w r and A * that the convergence of R α f is towards the projection P ϕ * of ϕ * onto N(A) ⊥ . In our case, the reconstruction calculates an approximation to P w * r . This completes the proof.
Usually, Tikhonov regularization is carried out by applying an efficient solver 1 to the equation which is equivalent to the spectral version of equation (39). Equation (42) is used for our numerical examples of the subsequent section.

Sensitivity Analysis
An important basic question is the influence of noise on the reconstruction. Here, we carry out a sensitivity analysis, i.e. we calculate the Fréchet derivative of the reconstructed kernel with respect to the input function u. Differentiability is obtained in a straightforward manner following [21], Chap. 2.6. We start with equation (35), where the operator V r and the right-hand side ψ r depend on the input function u. The reconstruction of w is carried out by the regularized version of which in the case of Tikhonov regularization (42) is We differentiate with respect to u on both sides and employ the chain rule and Eq. (2.6.21) of [21], to derive the unregularized form and the derivative of the regularized reconstruction where we use the notation The derivatives of V r and ψ r with respect to u are calculated as follows, where we restrict our presentation to the case where we are given one trajectory only. The operator V r in its dependence on u is given by leading to the Fréchet derivative where f denotes the derivative of the function f (s) with respect to its real argument s ∈ R. We need to assume that f is differentiable and that the derivative is continuous and bounded. The derivative of the adjoint V * r with respect to the L 2 scalar products on Ω and [0, T ], which is is given by for r ∈ Ω. We note that V * r is an operator into L 1 (Ω), which depends bounded continuously on r ∈ Ω. The Fréchet derivative of the function ψ r given by (26) is readily seen to be given by for (r, t) ∈ Ω × [0, T ]. We summarize the results in the following theorem. Proof Differentiability follows from the differentiability of all the operators in (46) following equations (46) to (52) of the above arguments.

Numerical Examples
The goal of this section is to demonstrate the feasibility of the inverse method for the reconstruction of spatial kernels based on the spatio-temporal neural field activity. We study the feasibility in Sect. 5.1 and the sensitivity with respect to variations in the input function u in Sect. 5.2.

Feasibility of Kernel Reconstructions
First, we consider a one-dimensional manifold embedded in a two-dimensional space, illustrating the method for a case with 10,000 degrees of freedom. Then an example involving a two-dimensional spatial domain evaluates the method for an inverse problem with more than 200,000 degrees of freedom for the kernel estimation. We first need to consider the role of boundaries in the neural field model equation (3) and its examples. For any distribution of neurons in space some activity u(r, t) depending on time t can be defined. Mutual influence in space is given by the integral in equation (3). In contrast to models based on partial differential equations, there is no direct boundary effect in these equations. However, • if one uses a local kernel w(r, r ) with strong connectivity only in a neighbourhood of r, boundary effects for neurons close to the boundary of the domain will appear, since less neurons are included in a neighbourhood there; whereas • if the activity of neurons close to the boundary is close to zero, usually such boundary effects remain negligible.
We will study a setup which avoids boundary effects by the choice of an embedding of a one-dimensional manifold into two-dimensional space in our first example, where there are always the same number of neurons in a neighbourhood of any neuron on the whole manifold. The second example instead limits boundary effects by using only small excitations close to the boundary in a two-dimensional neural patch.

Example 1
We start with a simple one-dimensional closed curve or manifold, respectively, embedded in a two-dimensional space. In particular, we study the dynamics of the activity field u(r, t) on the boundary ∂B R ⊂ R 2 of a disk with radius R, as displayed in Fig. 1. Here we consider that v = 1, and use a simple and smooth delay function for r = (x, y) and r = (x , y ) with r, r ∈ ∂B R based on the embedding into R 2 which is defined by This simple sandbox for testing our method hence can be considered as neurons growing on the boundary of a disk, but connecting directly through its interior. This is reminiscent of the thin exterior layer of grey matter containing neurons connecting through an interior bulk of white matter containing axons in the brain. However, we point out that this is a different setup from previous studies that superficially appear similar, where the spatial domain instead is a ring with periodic boundary conditions [38,39]. We implemented the delay neural field equation in MATLAB ® based on an Euler method for the time-evolution of the system with zeroth-order or first-order quadrature (rectangular rule or trapezoidal rule) for the integral parts of the integrodifferential equation. For the purposes of studying the kernel reconstruction on a Table 1 Parameter values for Example 1. Simulations have been carried out with N = 101 equally distributed nodes on the circle, N t = 50 time steps, and a time step size t = 0.2 for the inverse problem r 0 (cos(π ), sin(π )) σ 1.0 r 1 (cos(π/3), sin(π/3)) τ 1.0 For the one-dimensional example the kernel can be visualized as a two-dimensional scalar function w(r, r ). We display (a) the original and (b) the reconstructed kernel of the one-dimensional delay dynamics shown in Fig. 1 rather short temporal window this simple approach is completely sufficient and does not show any deficiencies compared to higher-order methods for the forward problem, as employed for example in [21,28,40]. We first solve the direct problem, i.e., calculate the time-evolution of the excitation field u(r, t). As initial condition, we choose the exponential function We prescribe a neural kernel of the form w r, r = c e −τ |r−r 1 | 2 e −τ |r −r 0 | 2 + e −τ |r−r 2 | 2 e −τ |r −r 1 | 2 for r, r ∈ ∂B R ⊂ R 2 with constants c > 0 and τ > 0. The full set of values used for our simulations are given in Table 1. This leads to delayed excitation of areas around three points r 0 , r 1 and r 2 equally distributed on a circle, where, with some delay, the excitation field around r 0 will excite the field around r 1 , the field around r 1 will excite the field around r 2 and the field around r 2 will excite the field around r 0 again. The function f is chosen to be sigmoidal as in equation (6). We have generated a classical oscillator, as can be seen in the snapshots in Fig. 1 (black curves). Its kernel is visualized in Fig. 2(a). Next we reconstruct the kernel by the inverse problem technique from the so obtained temporal evolution of the excitation field u(r, t) for some time window t ∈ [0, T ] according to equations (30) and (39). Given a discretized version of u(r, t) on nodes for = 0, . . . , N and k = 0, . . . , N t − 1, we calculate φ and ψ according to equations (25) and (26) and then employ the regularization (39) via (42) to solve equation (35) for r ∈ ∂B R . In Fig. 2, we compare the original with the reconstructed kernel in the case where no additional noise is added, carried out with α = 0.01 and find a very good agreement between both. As a test, we employ the reconstructed kernel with the same initial condition to calculate a reconstructed neural field u rec (r, t) on (r, t) ∈ ∂B R × [0, T ]. The original dynamics is shown in black in Fig. 1, based on the kernel (55) visualized in Fig. 2(a). The reconstructed dynamics is shown in red in Fig. 1, based on the reconstructed kernel visualized in Fig. 2(b). A very good agreement between original and reconstructed solution is observed.
Example 2 As a second example, we study oscillating two-dimensional neural field activity. Here, the dimension of the state space is higher with N = 21 × 22 = 462 spatial elements as shown in Figs. 3 and 4. Our approach is analogous to the onedimensional example, but now with 213,444 degrees of freedom for the possible connectivity values (see Table 2). We first simulate the neural field dynamics based on equation (3) Fig. 3. The kernel has been chosen to be of a form similar to equation (55), but now with points r 0 , r 1 and r 2 in the two-dimensional neural patch. This leads to an oscillating field in an area around these points r j with j = 0, 1, 2. The activation function f is chosen to be sigmoidal again. The initial condition is a Gaussian excitation around the point r 0 . For our simple tests, we again employ zeroth or first-order quadrature and Euler's method to carry out the simulation.
The kernel w(r, r ) with r, r ∈ Ω now lives on a subset U := Ω × Ω of a fourdimensional space, since Ω is a subset of a two-dimensional patch. Visualization of w(r, r ) can be carried out by either fixing r and showing a two-dimensional surface plot, or by re-ordering r and r into one-dimensional vectors, so that w(r, r ) can be displayed in full as a two-dimensional surface. The first approach is chosen in Fig. 4(c), where the white star indicates r . The second approach is shown in Fig. 4(a). Next, we solve the inverse delay neural problem and reconstruct the kernel based on equation (30) regularized as indicated by equations (39) and (42). Again, this is carried out by calculation of φ and ψ first according to equations (25) and (26), then solving equation (35) by regularization via equation (39) with the regularization parameter chosen as α = 0.1. This choice leads to a reasonable stability of the reconstructions combined with high reconstruction quality, and it has been chosen by trial and error.
Figures 4(c) and 4(d) display the original and the reconstructed kernel column, which represents the impact of the location at the black star to all other spatial locations of the neural patch. The result as displayed in (d) shows that the regular- ized reconstruction of the delay neural kernel is not perfect. However, it is working well if the field activity reaches specific parts of the neural environment. Otherwise the reconstruction is just zero due to missing input for the reconstruction equations and the regularization chosen here. The regularization penalizes the distance to the zero kernel function. Therefore, the results clearly demonstrate the feasibility of the method.

Sensitivity with Respect to Functional Input
In this section we will carry out a numerical sensitivity study of our first example to explore the dependence of the kernel reconstructions on the input function u. It complements our sensitivity analysis of Sect. 4.2.
We study the stability of the reconstruction when we add some random error to the measured signal u(r, t) displayed in Fig. 1. We remark that we need measurements of our signal which are differentiable with respect to time, since the calculation of ψ in (26) includes the temporal derivative of the signal. In practical situations, this would be achieved by a temporal smoothing of the signal. Here, for testing the sensitivity we have added a random shift of a temporally smooth signal in each of the analysis points. The amplitude of the signal is given by ε = 0.01, which corresponds to noise of 1% added to the measured temporal signal; compare Fig. 5. Now, we study reconstructions with different regularization parameters α, where larger α means we regularize in a stronger way, damping the error which comes from the measurement error. Figure 6 displays three different choices of α, where α = 1 leads to reasonable reconstructions, α = 0.1 shows kernel reconstruction still disturbed by noise, and α = 0.01 does not lead to satisfactory reconstructions at all.
According to Theorem 4.2 we have continuity of u → w, such that if we lower the error ε for fixed α, we need to have convergence to the reconstructed kernel in the case of no data error. Indeed, we obtain a figure similar to Fig. 6 when we lower the error parameter ε from ε = 0.01 to ε = 0.005 and ε = 0.001, leading to the reconstruction displayed in Fig. 2(b) for ε = 0. In the upper image we display the input signal u(r, t) independence of the point index of the discretized vector r and the temporal evolution t ∈ [0, T ]. The lower image shows the measurement error which has been added to the signal before a reconstruction has been carried out

Conclusions
The purpose of this work is to develop an integral equation approach for kernel reconstructions in delay neural field equations and to study its practical feasibility. We simulate the activity and evolution of a delayed neural field of Amari-type to develop an effective approach to reconstructing the neural connectivity. As a preparation for the inverse problem, this work includes an explicit study of the solvability of the direct problem of the delayed neural field equation (3). We provide an easily accessible functional analytic approach based on an integral equation and Banach's fixed-point theorem.
As our main result, we apply inverse problems techniques to reconstructing the neural kernel assuming that some measurements of the activity u(r, t) are given. We start by formulating a family of integral equations of the first kind. Since kernel reconstruction is ill posed, we need regularization to obtain stable solutions. As stabilization method we employ the Tikhonov regularization. A sensitivity analysis is carried out, showing that the mapping of the input u to the regularized kernel reconstruction is Fréchet differentiable. The derivative is explicitly calculated based on the integral equation approach.
Finally, we provide numerical examples in one-and two-dimensional spatial domains. These examples show that the regularized reconstruction of the delay neural kernel is practically feasible. We study the numerical sensitivity, by adding random noise of size ε (testing 1%, 0.1% and 0.01%) and studying the regularized reconstruction with different regularization parameters. In this work, we assume the delay function D to be given, as it would be the case when the delay is approximately proportional to the distance of the nodes under consideration. If D is unknown, w is known and u is measured, we can solve in equation (3) for u(r , t − D(r, r )) for all r, r and t. This is still ill posed, since it involves an integral equation of the first kind, but then the determination of D is reduced to the reconstruction of D from the knowledge of u(r , t − D(r, r )), which strongly depends on the form of the signal u and conditions we impose on D. If neither the delay D nor w would be given, the kernel K r of operator V r would be unknown and part of the reconstruction, leading to many open questions of feasibility and observability. In general, the reconstruction of both the kernel w and the delay D is an important nonlinear, far reaching and challenging problem of future research.
In summary, we have developed a stable and efficient approach for the reconstruction of the connectivity in neural systems based on delay neural field equations. We expect the approach to be extensible to a wide range of field models with delay, and in particular to be highly useful for analyses of experimental data in the domain of computational neuroscience. These methods allow for the reconstruction of the underlying 'synaptic footprint' of connectivity from available neural activity measurements, thus providing a basis for simulation and prediction of real phenomena in the neurosciences.