Integrals and series related to propagators of neural and haemodynamic waves

The propagator, or Green function, of a class of neural activity fields and of haemodynamic waves is evaluated exactly. The results enable a number of related integrals to be evaluated, along with series expansions of key results in terms of Bessel functions of the second kind. Connections to other related equations are also noted.


Introduction
Neural field theory (NFT) averages over the microscopic properties of neurons to obtain partial differential equations for mean levels of neural activity [1][2][3][4][5][6]. Propagation of activity across the cerebral cortex is governed by damped wave equations in NFT, whose Green functions, or propagators, embody the linear dynamics of such systems. In particular, the Green function represents the response evoked by a localized stimulus, which is a widely used method of probing brain function [7]. Such impulse responses have also been used to describe observed avalanche-like dynamics of neural activity [8,9].
Many useful cases of NFT propagators have a form that incorporates both damped wave propagation and regeneration of activity with a characteristic gain. Here, we analyse a particular form that applies to activity propagation across the twodimensional (2D) surface of the cerebral cortex, approximating it as flat. The three-dimensional integral that must be carried out to obtain the propagator has not been evaluated in NFT, although we discuss below how it can be related to propagators found in relativistic quantum theory. Its direct evaluation, as done here, puts the results in a suitable form for use in NFT and enables a variety of other integrals to be evaluated that do not appear to be included in standard tables such as [10][11][12][13][14]. Notably, the propagator that governs slow haemodynamic disturbances in brain tissue is a limiting form of the one for neural activity, although the physical mechanisms described are quite different. Hence, the present work also allows the haemodynamic response function to be evaluated exactly [15].
In §2, we outline the necessary NFT background and the form of the propagator to be evaluated. Section 3 is devoted to evaluating the relevant integral exactly, while §4 explores a variety of other integrals and series that derive from it. Finally, §5 provides a summary.

Neural field theory propagators
This section briefly provides the background to the central integral on which the subsequent material is based. This integral has arisen in various contexts in models of neural activity and resulting haemodynamics, so the same notation is kept here. Moreover, the descriptions are kept as close as possible to those in the previous work cited here to avoid any ambiguity or confusion.
Neural activity propagates as voltage spikes in axons, and NFT deals with spiking activity averaged over many neurons to form a field at scales of millimetres and up [2]. In a broad class of NFTs, the propagation of moderate-amplitude local mean neural spiking-activity perturbations ϕ across the 2D cortical surface from a delta-function stimulus at (r, t) = (0, 0) can be approximated as obeying the linear partial differential equation [1,2,6,7,16] 1 g 2 @ 2 @t 2 þ 2 g @ @t þ 1 À r 2 r 2 fðr, tÞ ¼ dðrÞdðtÞ, ð2:1Þ where r and t denote position and time, ρ is the characteristic range of the near-exponential distribution of axon lengths [1,[4][5][6] and γ is a characteristic damping rate. Equation (2.1) represents large-scale activity that propagates across the cortex as a damped wave at a speed v = γρ [1,6], and can be expressed in integral form via its propagator [1,2,4,5]. Equation (2.1) describes direct propagation of activity without regeneration by neural interactions. If we approximate the cortical surface as being flat and neglect boundary conditions, a Fourier transform of equation (2.1) implies the direct propagator [1] where k is the wavenumber and ω the angular frequency. This propagator is represented by the arrow at upper right of figure 1. If there is a mean gain g for regeneration of spikes at destination neurons, the outgoing activity at the next step is gL. Macroscopically, g is the differential change in outgoing ϕ per unit change in incoming ϕ at each stage of post-stimulus neural interaction, with g < 1 for stability [1]. The total activity T evoked by a delta input is the sum of local by the external stimulus, indicated by I in figure 1, and the terms on the right of figure 1 weighted by corresponding powers of g. This gives [16,17] Tðk, vÞ ¼ royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 211562 which converges for g < 1. Substitution of equation (2.2) into equation (2.4) gives where T 0 excludes local excitation (the n = 0 term in equation (2.3), which corresponds to the term I in figure 1) and has been divided by g for convenience; it represents the normalized propagating part of the response. From equation (2.6), it is evident that T 0 satisfies equation (2.1), except for the appearance of g, which leads to the unit term inside the square brackets on the left being replaced by 1 − g. Equations (2.3)-(2.5) exhibit the decomposition of neural propagation into the contributions of local activation (n = 0), direct propagation represented by L (n = 1), and propagation that involves intermediate interactions (n = 2, …). Therefore, T(k, ω) is the total propagator (or transfer function or Green function), and Lðk, vÞ is the direct propagator without regeneration [16], both of which depend on k only through its magnitude k.
If we take the limit g → 1, we obtain a wave equation that has been used to describe the cortical haemodynamic responses to localized stimuli [15]. These responses are responsible for the blood oxygen level dependent (BOLD) signal that underlies functional magnetic resonance imaging (fMRI). In this context, the Green function obtained from equation (2.1) without the unit term within the large parentheses at left is termed the haemodynamic response function [15].

Evaluation of the propagator
Because the transfer function in equation (2.6) only depends on k, T 0 is isotropic in 2D coordinate space. Fourier transformation of equation (2.6) then yields Using the property of the Dirac delta function that δ(ax) = δ(x)/|a|, and performing the u integral, we then find for r < vt, with T 0 (r, t) = 0 otherwise because waves cannot reach r > vt in time t.
We must now evaluate the series in equation (3.13), which has the form which has the same series expansion. Hence, for 0 ≤ g ≤ 1 and r ≤ vt, we obtain At large t, the dominant term on the right side of equation (3.18) varies as exp½Àgð1 À ffiffi ffi g p Þt, which approaches zero as t → ∞, given that we have assumed g < 1 for stability. Figure 2 shows examples of the propagator in equation (3.18) at different times for g = 0, 0.5 and 0.9. For g = 0, we see a decaying activity front propagating outward at v = γρ. At higher g, the early evolution is very similar, but as waves are regenerated, the small-r levels are higher and the overall decay is slower. These results improve on the approximate form recently obtained to illustrate T 0 (r, t) in [8], although the scaling arguments in that paper are not affected by this refinement.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 211562 It is worth noting that the modified version of equation (2.1) that includes g is a special case of the 2D telegrapher's equation [18]. If one writes ψ(r, t) = e γt ϕ(r, t), ψ obeys the Klein-Gordon equation for a scalar boson of mass ffiffi ffi g p , whose 2D solutions can be written in terms of Hankel functions [11,[19][20][21]. If the advanced ( past-propagating) part of the resulting propagator is discarded, the resulting retarded (future-propagating) propagator for ψ can be shown to yield equation (3.18) for ϕ. This approach, however, does not give the expressions for other integrals and series that we obtain via the present approach in this section and §4.

Related integrals and series
We can evaluate the ω integral in equation (3.2) by means of Cauchy's residue theorem, using a contour that follows the real axis and is closed clockwise by a semicircle in the lower half plane (the limit is taken in which the radius goes to infinity). This yields An extensive search of standard tables of integrals has not uncovered any expression for the integral in equation (4.1), but this result must be equal to the one in equation (3.18). Hence, we find royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 211562 for t > r/v and zero otherwise. The integral ð e Àivt K m ½a À ivb dv ¼ ffiffiffi ffi p p 2 m Gðm þ ð1=2ÞÞ a m b 2m e Àat=b ðt 2 À b 2 Þ mÀ1=2 , ð4:14Þ is then obtained by writing α = r/ρ and β = r/v so r = ρα, v = ρα/β and γ = α/β on the right of equation (4.13). An extensive search has not found either of equations (4.13) or (4.14) in standard tables. Taking the limit β → 0 in equation (4.14) yields a delta function on the left. The result can then be rearranged to give Although this result is somewhat baroque, such re-expressions of the delta function as limiting forms can be useful in derivations.

Summary and conclusion
The propagator (or Green function or impulse response function) of cortical large-scale neural activity and small-scale haemodynamic responses has been investigated, leading to the exact analytic form in equation (3.18). This expression reproduces several prior results in limiting cases, improves on a recent approximation, and enables the evaluation of a number of related integrals and series that do not appear to be contained in standard tables. In applications, the results will be of use in analytic work on responses evoked by stimuli, neural avalanches and analysis of fMRI scans. Relationships to retarded propagators of relativistic scalar bosons and the telegrapher's equation also exist but the direct derivation here is both instructive and yields a variety of new results.
Data accessibility. This article has no additional data.