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-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 two-dimensional (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-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]where r and t denote position and time, ρ is the characteristic range of the near-exponential distribution of axon lengths [1,4-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 . 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].

Figure 1

Diagrammatic representation of the terms in equation (2.3), illustrating that T is composed of local activation I plus terms representing direct propagation and indirect propagation via multiple neural interactions. A factor of g occurs at each vertex, representing regeneration of neural spikes, and the bare propagator carries activity between interactions.

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]which converges for g < 1. Substitution of equation (2.2) into equation (2.4) giveswhere T′ 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′ 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 (n = 1), and propagation that involves intermediate interactions (n = 2, …). Therefore, T(k, ω) is the total propagator (or transfer function or Green function), and 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′ is isotropic in 2D coordinate space. Fourier transformation of equation (2.6) then yieldswhere r = |r|. In all cases, causality requires that T(r, t) = 0 for t < 0, so we refer only to t ≥ 0 in what follows. Likewise, we assume that integrals extend over their whole domain unless otherwise indicated.

If we perform the k integral in equation (3.2), eqn 10.22.46 of [11] yieldswhere K0 is a Macdonald function; i.e. a modified Bessel function of the second kind. We next write the argument of K0 as λz withandsoThe multiplication theorem for Bessel functions of the second kind, or Macdonald functions, (eqn 10.44.1 of [11]) then implieswhich is valid for |λ2 − 1| < 1, a condition that is satisfied for g < 1.

Hence,The identity(eqn 10.32.8 of [11]) can then be used to rewrite equation (3.8) aswith v = γρ. Hence,becauseUsing the property of the Dirac delta function that δ(ax) = δ(x)/|a|, and performing the u integral, we then findfor r < vt, with T′(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 formand

Use of the identitythen yieldswhich has the same series expansion.

Hence, for 0 ≤ g ≤ 1 and r ≤ vt, we obtainAt large t, the dominant term on the right side of equation (3.18) varies as , 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′(r, t) in [8], although the scaling arguments in that paper are not affected by this refinement.

Figure 2

Propagator at various times for γ = 200 s−1, ρ = 0.05 m and v = γρ = 10 m s−1, showing T′(r, t) versus r from equation (3.18) for t = 5, 10, 15 ms, from left to right. (a) g = 0, (b) g = 0.5 and (c) g = 0.9.

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ϕ(r, t), ψ obeys the Klein–Gordon equation for a scalar boson of mass , whose 2D solutions can be written in terms of Hankel functions [11,19-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 yieldsAn 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

Some special cases of the propagator T′ are as follows:

To find the time-integrated activity that reaches a distance r, we integrate equation (3.2) over time, which is equivalent to extracting the ω = 0 component. Thus, [1,10,11,16]We can integrate equation (4.4) over r to obtain the total activity ultimately initiated by a unit impulse. This is equivalent to selecting the (k, ω) = (0, 0) component of T(k, ω), which is (1 − g)−1.

The result in equation (4.4) must be equivalent to the integral of equation (3.18) over t from r/v to infinity. Hence,or, equivalently,

If the total activity at a given time is sought, one must integrate equation (3.1) over position, which is equivalent to extracting the k = 0 component. This givesThis result can also be obtained by integrating equation (3.18) over r from the origin to r = vt. For g < 1, the right of equation (4.8) rises linearly before reaching a maximum and decaying as , whereas if g = 1 the initial linear rise is followed by an approach to the asymptotic constant level of . For g < 1, the integral of equation (4.8) over time is (1 − g)−1 which is always greater than 1 because of regeneration of activity.

If g = 0, we find from equation (3.18) which reproduces a result in [1]. This is also the expression for , as seen from equation (2.2).

A different family of integrals can be obtained by expanding the integrand in equation (3.2) in powers of g for g < 1. This giveswhere we have used eqn 2.12.4.28 of [13].

Since the value of g can be varied continuously, corresponding powers of g in equations (3.13) and (4.12) can be equated term by term. This yieldsfor t > r/v and zero otherwise. The integralis 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 giveAlthough 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.

Click here for additional data file.