Axonal Velocity Distributions in Neural Field Equations

By modelling the average activity of large neuronal populations, continuum mean field models (MFMs) have become an increasingly important theoretical tool for understanding the emergent activity of cortical tissue. In order to be computationally tractable, long-range propagation of activity in MFMs is often approximated with partial differential equations (PDEs). However, PDE approximations in current use correspond to underlying axonal velocity distributions incompatible with experimental measurements. In order to rectify this deficiency, we here introduce novel propagation PDEs that give rise to smooth unimodal distributions of axonal conduction velocities. We also argue that velocities estimated from fibre diameters in slice and from latency measurements, respectively, relate quite differently to such distributions, a significant point for any phenomenological description. Our PDEs are then successfully fit to fibre diameter data from human corpus callosum and rat subcortical white matter. This allows for the first time to simulate long-range conduction in the mammalian brain with realistic, convenient PDEs. Furthermore, the obtained results suggest that the propagation of activity in rat and human differs significantly beyond mere scaling. The dynamical consequences of our new formulation are investigated in the context of a well known neural field model. On the basis of Turing instability analyses, we conclude that pattern formation is more easily initiated using our more realistic propagator. By increasing characteristic conduction velocities, a smooth transition can occur from self-sustaining bulk oscillations to travelling waves of various wavelengths, which may influence axonal growth during development. Our analytic results are also corroborated numerically using simulations on a large spatial grid. Thus we provide here a comprehensive analysis of empirically constrained activity propagation in the context of MFMs, which will allow more realistic studies of mammalian brain activity in the future.

However, MFMs face severe technical difficulties when dealing with non-local neural activity, which is propagated across cortex by long-range axonal fibres. In order to incorporate the effects of such distributed activity a number of assumptions are typically made, the most important being a single value for the activity propagation delay between distant neural masses. This is the case even in otherwise sophisticated models, for example in those combining MFMs with Dynamic Causal Modelling (DCM) [24]. Most modelling approaches (e.g., [25,26]) follow here the lead of the seminal paper by Jirsa and Haken [27], who employed several simplifying assumptions to describe long-range activity propaga-tion with a partial differential equation (PDE). However, their ansatz still assumes a single value for the cortico-cortical axonal conduction velocity, and thus conduction delays between neural masses are exactly proportional to their distance with one uniform constant. We will show below that approximations made in deriving the actual propagation PDE result in an implicit velocity distribution, which nevertheless due to its origin remains strongly peaked at maximum conduction velocity and is onesided, i.e., there is an infinitely sharp cut-off at maximum speed. MFMs typically describe neural masses consisting of 10 5 to 10 7 neurons each. Thus even if the conduction velocity of one axon can be approximated well with a single conduction velocity, one should expect a distribution of conduction velocities between neural masses given the many axons involved. Empirical measurements of conduction velocities, either directly via conduction latencies or indirectly via fibre diameters, indeed suggest that conduction delays are rather broadly distributed. Initial attempts by Hutt and Atay [28,29] to incorporate broad axonal velocity distributions in a particular, spatially continuous MFM have revealed that such broad distributions maximize the speed of travelling front solutions. This may indicate the influence of natural selection optimizing information transmission in cortex.
Hutt and Atay [28,29] made use of a general integro-differential formula for activity propagation, which allows a straightforward introduction of velocity distributions. It is just this integrodifferential formula, which is commonly simplified towards a PDE [27]. As discussed for example by Liley et al. [26], local PDE formulations offer a number of significant advantages over their non-local (integral) counterparts. In particular, they enable the use of powerful analytical and numerical analysis methods, at least for specific spatial wavenumbers, and allow the application of standard numerical techniques for the solution of MFMs. The latter point is particularly important for large-scale simulations, see for example [9,13], where computation speed is essential. As derived in [30] by the present authors, one can always extract the velocity distribution implied by the PDE formulation of an MFM. But so far the exact form of these distributions have been largely an accidental side product of approximations. It is hence no surprise that the velocity distributions of models in current use are unsatisfactory. Incorporating a sensible velocity distribution into an analytically and numerically tractable PDE formulation has not been achieved before.
Motivated by physiological and anatomical fidelity on one hand, and by computational necessity on the other, we here introduce a novel PDE formulation describing the propagation of cortico-cortical axonal activity that incorporates monotonically decaying synaptic connectivity with a smooth unimodal distribution of axonal conduction velocities. We obtain good fits with our new model to experimental data on conduction velocities derived from myelinated fibre diameter measurements in the human corpus callosum [31]. This allows for the first time to simulate long-range conduction in humans based directly on experimental findings. A straightforward extension of initial propagator ansatz also allows us to fit data from lower mammals, which generally feature less small diameter (myelinated) fibres. Studying activity conduction in animal cortex is important in its own right, but also significant for the suitability of animal models for human studies. For example, the CoCoMac database [32,33] contains precise information on the connectivity of macaque cortex from extensive tracer studies, which cannot be obtained similarly from humans since such techniques are lethal. While CoCoMac connectivity can be mapped to human cortex [34] and calibrated with human connectivity data from non-invasive Magnetic Resonance Imaging [35], the question would remain whether similar anatomical connections actually serve the same function. Clearly an improved understanding of the dynamics of activity conduction in animals and humans is of great significance to this question.
We obtain reasonable fits with our extended ansatz to extensive unmyelinated and myelinated data from rat subcortical white matter [36], and discuss briefly the clear differences that exist to the fit to human callosal data. Finally, we also analyse analytically and numerically the dynamical impact of using our new propagator. Following the methods in Coombes et al. [30], we can show that in contrast to the most commonly used longwavelength propagator, our realistic velocity distributions enable the formation of spatio-temporal patterns for smaller perturbations in mean neuronal firing rates. This may follow more closely the biological situation, where a range of energetic constraints need to be negotiated in order to ensure that pattern formation, and thus perception, occurs in metabolically optimal circumstances. We confirm these results with some explorative computational simulations on large spatial grids using our novel propagator. So far, conduction parameters in mean field models have been either chosen largely arbitrarily from a wide range of plausible values, or adjusted freely to help reproducing the phenomena under investigation. Our fits to human and rat data, and future fits to other experimental data using our methods, constrain propagation parameters empirically and independently. This will reduce considerably the uncertainties of future predictions using the mean field framework.

Dispersive propagator
In most neural field models developed to date the activity variables that are spatially propagated are the local mean neuronal population firing rates, S j . Because action potentials propagate with a finite conduction velocity, the mean rate of arrival of presynaptic impulses w jk to cells of type k from neurons of type j can be written as a time-retarded integral of the respective distant local mean excitatory neuronal firing rates: ð ?
where spatial integration occurs over a two-dimensional planar cortical sheet C (x,x 0 [R 2 ). The distance-dependent velocity distribution function f jk v jx,x 0 ð Þtakes into account that fibre paths with different conduction velocities can exist between different domains. This conditional distribution is normalised such that Ð ? 0 dvf jk v j x,x 0 ð Þ~1. The function w jk x,x 0 ð Þ is the synaptic footprint that describes the geometry of network connections. The distance dependent Green's function, G jk , is defined as: In the absence of detailed anatomical data it is common practice to consider synaptic connectivity functions to be homogeneous and

Author Summary
Due to the sheer number of neurons and the complexity of their interactions, the modelling of brain activity is particularly challenging. How can computationally tractable models of brain function be developed that are nevertheless biologically plausible? The ''mean field'' approach, borrowed from statistical physics, is to model the average activity of populations of neurons rather than the behaviour of individual neurons. While a large number of promising theories have been developed with this approach, they fall short of biological fidelity in the way interactions between distant populations have been modelled. In particular, it is often assumed that all neurons interact via connections of very similar conduction velocity, when in fact experiment suggests quite the opposite: populations of neurons are connected by axonal fibres with a broad range of velocities. We develop here activity propagators that provide for the first time the ability to realistically and efficiently simulate connectivity in mean field theories, and demonstrate how to use them to fit successfully experimental data from both human and rat. With our novel propagators, one can thus study on an empirical basis the role of activity propagation in both healthy and diseased mammalian brains.
isotropic so that w jk x,x 0 ð Þ:w jk x{x 0 j j ð Þ . We will also assume that this restriction applies to the velocity distribution functions, i.e., f jk v jx,x 0 ð Þ:f jk v | x{x 0 j j ð Þ , and therefore G jk x,x 0 ,t{t' ð Þ : G jk x{x 0 j j,t{t' ð Þ . This assumption of isotropy can be relaxed at the price of increased computational effort [30,37,38], as will be discussed below in a separate subsection. The right hand side of (2) now has a convolution structure, and its Fourier transform, where k~k j j. If G jk (k,v) has the form R jk (k 2 ,iv)=P jk (k 2 ,iv) then the integro-differential Eq. (2) can be written as the equivalent PDE R jk ({+ 2 ,L=Lt)w jk x,t ð Þ~P jk ({+ 2 ,L=Lt)S j x,t ð Þ, i.e., the corresponding partial differential operators are obtained with the Fourier replacements k 2 ?{+ 2 and iv?L=Lt.
The most common propagator form used in mean field models of electroencephalographic activity derives from the following simple ansatz for the Green's function: an exponential decay with distance of propagated firing rates is combined with isotropic conduction where r:jxj §0 and axonal velocityv v jk w0 together imply the causal conduction of activity through a Dirac d distribution of delays. The normalization constant w 0 jk counts the total number of synaptic connections made by the axonal fibres originating from neurons of type j that terminate on neurons of type k. The exponential decay with the characteristic distance scaleŝ s jk should be understood as due to diminishing connectivity [39], rather than as decay of the amplitudes of the action potentials themselves. The Fourier domain propagator in Eq. (5) is non-polynomial, but can be approximated for small k, and hence long wavelengths l~2p=k, with a polynomial form. Settingṽ v jk : ffiffiffiffiffiffiffi ffi 3=2 pv v jk ands s jk : ffiffiffiffiffiffiffi ffi 3=2 pŝ s jk we obtainG where H is the Heaviside step function, which now maintains causality. We will subsequently refer to this as the long-wavelength approximation. The standard inhomogeneous, 2-dim. telegraph equation [25][26][27] results Note that (7) is a special case. If we substitute then jk obeys an inhomogeneous wave equation. Note that Eq. (8) corrects a sign error in Eq. (61) of Ref. [25]. The approximate impulse responseG G jk (r,t) in Eq. (6) can hence be recognized as that of a 2-dim. wave with velocityṽ v jk multiplied by an exponential decay with velocity-dependent distanceṽ v jk t.
The infinitely precise conduction delay d t{r=v v jk À Á of ansatz Eq. (5) is at odds with the broadly distributed delays measured by experiment. In the next section we will show that the longwavelength approximation largely inherits this problem. An obvious amelioration would be to use a Gaussian normal distribution of delays: where c is an appropriate normalization constant and the Heaviside H enforces causality. However, Eq. (9) leads to the same type of nonpolynomial Fourier structure as Eq. (5), only multiplied with exp½(ivs delay ) 2 =2. Thus again an approximation would be needed to obtain a polynomial form and hence a PDE. A key observation is that the problematic fractional power 3=2 arises from the spatial Fourier transform of exp ({ar) terms, where the a are independent of distance but can depend on time, and that we can eliminate all such terms from the ansatz by setting s delay ? ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ts jk =v jk p : We can Fourier transform this expression, first spatially (which is equivalent to a zeroth order Hankel transform) and then temporally, even if it is multiplied with powers of t. Hence we now propose the following Green's function: where nw0 and C(n) is the Gamma function with C(n)~(n{1)! for integer n. The corresponding Fourier domain propagator is Using this to propagate local mean firing rates according to Eq. (4) is hence equivalent to the following two-dimensional PDE where only n[N 1 will realize any practical benefits for analysis and computation. Note that for n~1 this corresponds to a twodimensional, inhomogeneous cable equation. We will subsequently refer to this novel ansatz as the dispersive propagator. It should be emphasized at this point that single propagation PDEs, like the dispersive Eq. (13) and the long-wavelength Eq. (7), imply that firing rate activity passes continuously between any two arbitrarily chosen cortical locations. However, cortico-cortical fibres are known to also selectively connect separated areas of cortex in a direct manner, see for example Ref. [40]. Such nonlocal propagation cannot be modelled with the PDE descriptions of activity conduction described so far. To include non-local effects one must either resort again to the general integral equations, or map cortex to a mixture of overlapping patches based on a chosen PDE description. Recently good progress has been achieved for the latter option [38], in particular also by turning such descriptions into a kind of DCM [41], which makes possible robust fits to experimental neuroimaging data. Our efforts here are complementary to these pioneering works, since we are concerned with obtaining physiological conduction velocity distributions in the typical PDE framework. For example, the long-wavelength approximation Eq. (5) in Ref. [38] could be replaced with our dispersive Eq. (13) as basis for considering non-local effects, thereby increasing the realism of the non-local conduction model even further. We will explain in a separate subsection below in what way anisotropy and inhomogeneity can also affect the extraction of velocity distributions from experimental data.
In the original ansatz of Eq. (5), impulses would arrive at distance r from a source precisely after a time r=v v jk had passed. The extension in Eq. (9) was constructed such that the impulses would arrive with a Gaussian normal distribution of delays having mean r=v jk and standard deviation s delay . We can recover this from the respective Green's functions by computing the statistical characteristics of delays, appropriately normed by the decay of connectivity to distance r: StT: , Thus indeed StT~r=v v jk and s t~0 for the original ansatz Eq. (5), but for the long-wavelength approximation Eq. (6) thereof one finds instead From the results for r=StT one can see that the characteristic longwavelength (ṽ v jk ) and dispersive (v jk ) velocities still indicate the axonal conduction velocities, but only on average and at large distances. A ''large'' distance means here one much greater than the characteristic decay scales of connectivity,s s jk and s jk , respectively. At large distances the standard deviation of delays s t becomes constant for the long-wavelength approximation, but s t * ffiffi r p for the dispersive propagator, i.e., it grows with the square root of distance. We also see that s t ?
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi StTs jk =v jk p at large distances, which recovers the substitution of s delay leading to Eq. (10). Note finally that long-wavelength StT is identical to the dispersive StT at all distances forṽ v jk~vjk and n~1.

Synaptic connectivity and velocity distribution
By integrating the dispersive Green's function Eq. (11) over time we obtain the implied dependency of synaptic connectivity with distance w 0 Here w 0 jk counts the total number of synapses formed and w jk (r)=w 0 jk is the probability distribution of the synaptic footprint, i.e., the likelihood that a synapse forms at distance r, where Ð ? 0 dr2prw jk (r)=w 0 jk~1 . Connectivity w jk (r) remains finite for r?0 only if nw1, in which case w jk (0)~w 0 jk =½4p(n{1)s 2 jk . In practice the nƒ1 divergence for r?0 is of little concern, as neural field models are not meaningful below some minimal size r avg over which mean population activity is defined. The contributions of synaptic connections within the disc 0ƒrƒr avg to the total number of synaptic connections w 0 jk vanishes for r avg %s jk for all nw0. Eq. (19) should be compared with the connectivity function for the long-wavelength approximation Eq. (6) withw w jk (r) normed to w 0 jk as in Eq. (20). We note again an equivalence to Eq. (19) with n~1.
Both the dispersive and the long-wavelength propagator thus have synaptic footprints decaying with distance *r n{1 K n{1 r=s jk À Á , where for the latter n~1. However, experimental counts of synaptic connectivity usually have been fit with the simpler exponential decaŷ Thus the question arises whether dispersive connectivity is compatible with data that apparently fit an exponential decay, and whether one can use such previous fit results to constrain also the dispersive propagator. An exponential decay is also what the original ansatz Eq. (5) used. Hence in previous works it has been assumed that model and fit scale are basically the same quantity. But it will become clear now that after the long-wavelength approximation Eq. (6) this is not correct anymore. Let us assume that the dispersive synaptic footprint Eq. (19) with parameters w 0 jk and s jk represents the true underlying distribution of connectivity, and that from it parametersŵ w 0 jk :dw 0 jk andŝ s jk :es jk are estimated with a fit assuming the exponential distribution Eq. (22). Therefore we wish to determine which d and e best corrects for the mismatch. In practice, experimental counts of synaptic connections are usually sorted into distance bins ½r i ,r iz1 , where r i~i : Dr with i~0, . . . ,i max . We can scale r=ŝ s jk :x and r=s jk :ex, whereŝ s jk is known from the experimental fit. The counts per bin are then and we can minimize this expression explicitly to determine d and e. To give a numerical example: assume i max z1~20 bins of width Dx~x iz1 {x i~0 :25, i.e., the bin size was a fourth of the fittedŝ s jk and in the last bin connectivity had decayed to less than one percent of maximum. For different powers n we can then obtain numerically scaling factors d and e: n d We find that the normalization correction d has an asymptotic value for large powers n, whereas the decay correction e grows as ffiffi ffi n p . The resulting synaptic connectivity is shown in Fig. 1A. For simplicity we have assumed here thatŵ w 0 jk~2 pŝ s 2 jk , i.e., that w w jk (r)~exp ({r=ŝ s jk ). The dispersive curves are hence w jk (r~es jk x)~e 2 =½2 n{1 C(n)d : (ex) n{1 K n{1 (ex) with the scaling factors derived above. While we show continuous curves here, the correction was performed for binned data. It is obvious from the reasonably close match that dispersive connectivity may well be mistaken for an exponential decay, given the large statistical and systematic errors typically involved in synaptic counts. Note that the n~1 divergence for small distances would not be visible in a binned count. Nevertheless, it is obvious that the n~1 case, and hence the long-wavelength approximation, does not match an exponential decay better than higher powers of n. Furthermore, ð26Þ for n~1 in this example we find the optimal scaling w 0 jk~1 :673 :ŵ w 0 jk and s jk~2 :385 :ŝ s jk . In general for long-wavelength models one should actually choose w 0 jk and s jk which are significantly larger than those measured in experiments. Note that our long-wavelength decay scale absorbed an expansion factor ffiffiffiffiffiffiffi ffi 3=2 p to keep Eq. (6) simple. Without this, scaling by ffiffiffiffiffiffiffi ffi 2=3 p : 2:385^1:947 would be required here.
Eq. (3) enables us to determine the underlying conduction velocity distribution of the axonal fibres that arises from our newly proposed dispersive propagator. Thus we obtain Using the Green's function Eq. (11), the distance-dependent velocity distribution f jk v jr ð Þ becomes which has a maximum at The f jk v j r ð Þ distribution indicates the probability of conduction velocity v at a given distance r. As far as experimental data are concerned, this distribution is appropriate for measurements of conduction latencies between brain regions. For that case we can consider r to be fixed and note that f jk v jr ð Þ is properly normed as a conditional probability distribution in v, i.e., Ð 1 0 dvf jk v j r ð Þ~1. The time t~r=v max indicates the moment when most propagated activity arrives at once in a region. One can speculate that this has the highest likelihood to induce a signal visible over local background activity. According to the first limit in Eq. (29), we then expect latency data for distant (r&ns jk ) regions to measure conduction velocities *v jk . Fig. 2A shows a plot of the cumulative distribution corresponding to Eq. (28). We prefer to show the cumulative distribution here, because of the large variations of f jk v j r ð Þ in the shown range of v and r. Furthermore, this allows a direct comparison with the long-wavelength approximation later on. The sigmoidal shape of F jk v j r ð Þ in v corresponds to the unimodal form of f jk v j r ð Þ. The position of the mode v max of f jk v jr ð Þ is indicated by a solid black line on the F jk v j r ð Þ surface. That F jk v max j r ð Þv0:5 indicates that the distribution is skewed towards higher velocities. However, we can see that the distribution becomes less skewed for larger r. Furthermore, we see that neuronal populations at greater distances on the cortical surface are connected by faster fibres. While from a functional perspective this makes intuitive sense, there is at present no direct anatomical or histological evidence for this. We discuss some indirect evidence below. The second limit in Eq. (29) shows that higher order n distributions describe overall slower connectivity for the same v jk .
The distance-dependent connectivity function for each fibre system of velocity v, w jk (v,r), is then where Ð dr 2pr Ð dvw jk (v,r)~w 0 jk , and w 0 jk counts the total number of synapses formed. Hence, w jk (v,r)=w 0 jk defines the joint probability distribution for propagation with speed v to distance r. The marginal propagation velocity distribution over all r is then where Ð ? 0 dvf jk (v)~1. As far as experimental data are concerned, this distribution is appropriate for measurements of local fibre diameter statistics, which can be related to conduction velocities. Such statistics catalogue all fibres passing through a local slice, irrespective of the distance between the neural populations they connect. This corresponds to integrating out the distance in Eq. (32).
We show the marginal velocity distribution (multiplied by the constant v jk ) in Fig. 1B for several different powers n. The rapid sharpening up of the distribution for higher powers is readily apparent. The statistical characteristics of the dispersive f (v) distribution are collected in Tab. 1; note also that it becomes a beta-prime distribution with a~1 and b~n under nonlinear scaling x:v 2 =v 2 jk . For nƒ1=2 both the mean and standard deviation of the dispersive f (v) do not exist, like for a Cauchy random variable, and for 1=2vnƒ1 the mean exists but not the standard deviation, due to the tail-thickness of the distribution. Thus at n~1 large variations of the conduction velocity are probable. The coefficient of variation s v =SvT asymptotes to 0:523, even then indicating a broad distribution. For n~2,3,4 the corresponding velocity distributions already have 66%, 79% and 84%, respectively, of this maximal ''sharpness''. Skew c 1,v exists for nw3=2 and indicates preference for higher velocities. The mode v mode of the marginal dispersive velocity distribution is smaller than v jk , see Tab. 1. This is more pronounced for higher order n due to a larger fraction of slower fibres. By contrast, the mode v max of the conditional dispersive velocity distribution approaches v jk for large distances, see Eq. (29), but again more slowly for larger n. Both mode speeds are identical in the dispersive case for r~s jk ffiffiffiffiffiffiffiffiffiffiffiffi 1z2n p , where below this distance v max vv mode and above this distance v max wv mode . As we see from this example, comparisons of the dominant speeds -v mode estimated from fibre diameters in a local slice and v max from latencies between distant brain regions -can be used as an experimental probe of the underlying connectivity. For fibre distributions like the dispersive one, in which more distant regions are connected by faster fibres, one would expect distance-dependent relations between v mode and v max qualitatively similar to the ones just described. Latencies observed at different distances could complement the experimental constraints from local fibre diameter measurements quantitatively, too. However, the difference between v mode and v max becomes more significant for measurements at larger distances, where unfortunately one would also generally expect worse signal to noise ratios. Thus it is currently unclear whether such comparisons are in fact feasible experimentally beyond a qualitative consistency check. Nevertheless, there is a chance to gain significant new insights into brain connectivity here using comparatively ''simple'' techniques, or even from a re-analysis of previously obtained data.
The distance-dependent velocity distribution for the longwavelength approximation Eq. (6), unlike for the dispersive propagator, is truncated for velocities greater than v jk : Again for r fixedf f jk v jr ð Þ becomes a conditional probability distribution in v appropriate for comparisons with experimental conduction latencies. Fig. 2B shows a plot of the corresponding cumulative distributionF F jk v jr ð Þ, integrated as in Eq. (30). Note thatf f jk v jr ð Þ?? for v?ṽ v jk , whereasF F jk v j r ð Þ is well-behaved in the limit and hence can be plotted easily. For rv For very small r this maximum even formally becomes dominant, but at such distances the MFM loses validity. Thus the global maximum is in practice always determined by the cut-off v v max~ṽ v jk . The position of the maxima off f jk v j r ð Þ is indicated in Fig. 2B by two solid black lines on the surface ofF F jk v jr ð Þ. The corresponding marginal velocity distribution, which can be related to measurements of axonal diameters, is given bỹ and its statistical characteristics also are collected in Tab. 1. We see that this distribution is very sharp, with a coefficient of variation s v =SvT~0:284, and skewed to lower velocities. Indeed, high velocities are cut off atṽ v jk . Note that the mode of the marginal distribution is the sameṽ v jk as the maximum velocity between distant brain regions of the distance-dependent distribution. Thus here we would predict that fibre diameter and latency measurements derive roughly the same conduction velocity. Basically the long-wavelength approximation retains the original sharply peaked velocity distribution of fibres with a single conduction velocityṽ v jk . If the comparison between conduction velocities derived from diameter measurements and latencies can achieve sufficient statistical significance, then this would allow an experimental distinction between the dispersive and long-wave-length propagators. We consider investigating inter-hemispheric connectivity between contra-lateral brain regions as promising, because it is heavily dominated by just one fibre type (myelinated fibres), with fairly homogeneous regional expression across large distances. This adds particular significance to our fit of diameter data of myelinated axons from human corpus callosum performed below.
Incorporating anisotropy and inhomogeneity. In our presentation of the dispersive propagator, and the subsequent derivation of the conditional and marginal velocity distributions, we have assumed both isotropy and homogeneity of the corresponding connectivities. It is fortunate that these restrictions can be relaxed, given that neither homogeneity nor isotropy would be expected to hold fully in real brains, particularly not so for long-range connectivity. First, inhomogeneities will be described well by our equations in an average sense, as long as they are relatively small and random according to some unimodal distribution, e.g., a normal distribution. This fits well with the general MFM approach of describing only the ''mean fields'' of cortex. Further, the parameters may vary in an arbitrary inhomogeneous fashion over distances farther away than a few times the characteristic scale of synaptic connectivity s jk , without causing local complications. Conducted over these distances a local pulse will have mostly decayed away, hence the PDEs remain valid. This suggests a separation of cortex into regions of ''homogeneous enough'' conduction properties. If the inhomogeneous variation of conduction properties across cortex is nevertheless smooth, then even a single PDE with matching spatial variation of parameters Table 1. Statistical characteristics of the dispersive, long-wavelength, and difference marginal velocity distributions.
:8660 Statistics are shown for the following marginal velocity distributions: dispersive Eq. (32), long-wavelength Eq. (36), and difference Eq. (51). The characteristic velocities v v for these three distributions are v jk ,ṽ v jk , and v 1 , respectively. SvT, s v , and c 1,v are the mean, standard deviation, and skewness in v, respectively. In order to achieve a compact notation, we have defined g(n):C n{ 1 , and l c : . For w 2~0 one finds w r~1 , as the difference propagator turns into the dispersive one. We have not found a closed analytic form for the mode v mode and median v median of the difference distribution, but they can be computed numerically. Further definitions needed for the evaluation of the difference distribution statistics are collected in Eqns. (43) and ( could be used as model. Otherwise one would have to take special care at the boundaries. Second, to describe anisotropic conduction a generalization to ''patchy'' propagators is possible. Work by Robinson [37] has shown that one can generate basically arbitrary angular modifications of conduction properties at the price of introducing more PDEs. Basically this technique relies on a spatial Fourier decomposition of long range connectivity. Hence the sharper the anisotropy one wishes to describe, the more PDEs one has to employ. See for example Ref. [30], where sinusoidal variations in two principal directions required the solution of four coupled complex PDEs, instead of one real PDE. In practice a compromise between biological fidelity and numerical complexity has to be made. Consider then the following ''patchy'' Green's function which is homogeneous but anisotropic. It allows the specification of anisotropic connectivity through a decomposition into an isotropic Green's function G jk and an anisotropic, but timeindependent, modifier M jk . Now we can use Eq. (3) for G M jk and integrate over the Dirac d-distribution, as for Eqns. (27) and (28), but without any assumption of isotropy. The synaptic footprint is again the integration over time of G M jk , like in Eq. (19). Thus the conditional velocity distribution becomes here i.e., the anisotropic modifier M jk cancels out and the conditional velocity distribution f M jk is found to be isotropic, and identical with the f jk of the isotropic Green's function G jk . Thus an isotropic conditional velocity distribution is entirely compatible with anisotropic connectivity.
Rewriting Eq. (3) in polar coordinates, x{x 0 j jand h, one finds that in general and thus the potential anisotropy of propagation velocities is independent of any evidence or assumptions regarding the anisotropy of synaptic connectivity. In other words, how fast the fibres connecting two regions are is a different question to the number of fibres that connect these two regions. Hence even for real brains one can start with the parsimonious isotropic assumption for the conditional velocity distribution Þ , and assume that anisotropies are due only to w jk h, x{x 0 j j ð Þ . Then the fibre system is potentially anisotropic, but where fibres grow their distribution of conduction velocities is not dependent on the direction in which they are growing. Further, define the ''angular average'' Then the generalization of Eq. (32), which assumes that the conditional velocity distribution is isotropic but allows for anisotropy in the connectivity, can be written as where we have set r~jxj again. This clearly depends only on w ShT jk (r), and may very well be practically indistinguishable from isotropic conditions. For example, a fibre system with one strongly dominant direction w jk (h,r)~w jk (r)d(h{h 0 ), which is roughly the case within corpus callosum, yields the same isotropic f jk (v) through the renormalization of w 0 jk . For these reasons we will continue with the assumption of isotropy for fits of the marginal velocity distributions to data. However, more precise data on both connectivity and conduction latencies may well make possible in future to disentangle anisotropies further, potentially showing that our parsimonious assumption of an isotropic conditional velocity distribution was incorrect. One also needs to keep in mind that for simulations of cortex the introduction of inhomogeneous regions and ''patchy'' propagators will likely be required to achieve good biological fidelity, even if one assumes isotropic velocity distributions. In this regard the methods of Daunizeau et al. [38] may prove particularly useful, which systematically map conduction PDEs to heterogeneous cortico-cortical connectivity in the human brain.

Difference propagator
Finally, there appears to be a general trend in experimental data that higher mammals have a larger proportion of small diameter fibres, see for example the discussion in the section ''Species differences'' of [31]. We will encounter this phenomenon when trying to fit human [31] and rat data [36]. Small fibre diameters correspond to low conduction velocities, as we will see in detail below. Unfortunately the dispersive propagator predicts too much low velocity conduction, and thus a too large fraction of small diameter fibres, to fit the rat data well. Whereas the longwavelength approximation fails entirely to describe either human or rat data, but because of high, not low, velocity conduction: its marginal velocity distribution is sharply peaked close to an upper velocity limit, while all data require a broad, unimodal velocity distribution. We have been unable to find another single propagator equation, which both yields the polynomial Fourier structure leading to a PDE and describes the data from lower mammals better.
A constructive approach for dealing with this problem posed by animal data has however proven successful. The basic idea is to subtract two dispersive propagators w tot~w1 {w 2 , where the second dispersive propagator conducts activity more slowly, so that the resulting distribution is reduced at small velocities. This construction we will then call the difference propagator. Before we provide further mathematical details, we wish to justify this method with regards to the actual biology it is supposed to describe. Clearly there are no ''anti-fibres'' in the brain, hence w 2 and therefore also w 1 lack any direct biological meaning taken separately. But the biological meaning of the constructive solution w tot is not necessarily compromised, since in the end it is actually w tot which is compared with empirical measurements. The dispersive and long-wavelength propagators we have investigated so far are biologically meaningful and appropriate because of the following characteristics: First, they correspond to a Green's function non-negative for all positive times and distances. This implies that a positive local pulse also leads to positive pulses arriving at distant synaptic terminals. The impact of these pulses may be ''negative'', if they excite inhibitory populations, but the action potentials themselves do not somehow change sign. Second, synaptic connectivity has a roughly exponential decay with distance, as is appropriate for describing background connectivity in the brain. Third, the distance-dependent velocity distribution has a dominant mode, i.e., there is a preferred conduction velocity leading to typical latencies between brain regions. Fourth, the marginal velocity distribution has a shape which compares favourably with fibre diameter distributions. We will construct our difference propagator so that it shows all these characteristics. Hence while it may be less intuitive, and requires more computational effort, w tot will be as valid in terms of biology as the dispersive and long-wavelength propagators.
We first compute the ratio of two dispersive Green's functions G jk (r,t) from Eq. (11), which have different parameters : with normed spatial variables t:v 1 t=(2s 1 ) and y:r=(2s 1 ). The inequality is valid for powers n 1 ,m[N 1 , and thus n 2 wn 1 , as well as factor 0vf ƒ1, and we have set z:e m n n1 If we now define then it is clear that for w tot local firing S will be propagated with a combined Green's function G tot (r,t)~G 1 (r,t){G 2 (r,t) §0. By construction we have made certain that no unbiological ''negative pulses'' can arise here in spite of the subtraction. Thanks to the linear combination, the distributions are computed trivially, e.g., synaptic connectivity is Note that as integral over G tot , see Eq. (19), w tot (r) and hence w 0 tot must be positive, since G tot §0 and not zero in the entire integration range. In practice w 0 tot is the biological quantity and determines w 1 via Eqns. (48) and (44). We can again compute how this connectivity compares to an assumed exponential decay, as explained at Eq. (25). The sum to be minimized becomes now where r=s 1 :ex. For different powers n we obtain here scaling factors d and e which are similar to those of the dispersive propagator: n In Fig. 3A the corresponding difference connectivity is shown. We see that it may become feasible to measure experimentally the deviation to an exponential decay in particular for small x and high powers n, though overall the shape is still roughly exponential. The distance-dependent velocity distribution f tot v j r ð Þ~rG tot (r,t~r=v)=½v 2 w tot (r) and the distance-dependent connectivity w tot (v,r)~rG tot (r,t~r=v)=v 2 are of course also positive. It is straightforward to show that for r=s 1 &1 the conditional distribution f tot v j r ð Þ is indeed unimodal, with the maximum given by Eq. (29) upon replacing n?n 1 , s jk ?s 1 , and v jk ?v 1 . At r~2s 1 the mode velocities of the dispersive and difference propagators already differ by less than 10% for n 1 w1.
Finally we can compute the marginal velocity distribution Its statistical characteristics can again be found in Tab. 1. As before mean SvT only exists for n 1 w1=2, standard deviation s v only for n 1 w1, and skewness c 1,v only for n 1 w3=2. No further condition is required, since n 2 wn 1 . We have not been able to find analytic expressions for v mode and v median for unspecified powers n 1 and n 2 . However, computing them numerically for chosen powers is straightforward. Since we wish to deplete f tot (v) at small v, we want to maximize positive skewness c 1,v using the still available factor f . There is a clear mode of c 1,v in the range 0vf ƒ1, but again it is too difficult to find it analytically. Instead we obtain 225 numerical solutions for n 1~2 , . . . ,16 and m~n 2 {n 1~1 , . . . ,15, and then obtain a good three parameter fit for maximum skewness: With Eq. (52) we complete the specification of our difference propagator. In practice then, the difference propagator can be computed using two PDEs where only four parameters are actually free: n 1 , s 1 , v 1 , and n 2 . All the other parameters are dependent, see Eqns. (43), (44) and (52). Furthermore, n 2 wn 1 is also required. Thus in comparison to the dispersive propagator only one additional parameter is introduced here: the chosen power n 2 of the subtracted dispersive propagator. While the variables w 1 and w 2 have no independent meaning here as such, both describe independently propagated activity since their PDEs are not coupled. Hence one can think of w 1 as representing a ''full'' propagator, which one would encounter in humans, and of w 2 as representing a ''depletion'' propagator, which then removes activity conduction lacking in lower mammals.
Comparing dispersive and difference distributions for n 1~n in Tab. 1, we find now that both mean and standard deviation of the difference distribution are larger, but its coefficient of variation is smaller. Thus the difference distribution is sharper. Skewness is indeed more positive for the difference distribution, indicating the increased preference for higher velocities we aimed for. For m?? one finds w 2 ?0, i.e., the difference distribution becomes the dispersive one again. The m~1 case then also turns out to be least similar to the dispersive one concerning statistical characteristics. Our skewness fit cannot be expected to be faithful outside of the fit range, which however is sufficient for all practical purposes. The only exception is n 1~1 where skewness does not exist, but which may be of interest. The approximation in Eq. (52) extrapolates viably in that case with 0vf ƒ1, and for simplicity's sake we adopt here the fit for all n 1 . The resulting distribution is shown in Fig. 3B. In comparison to Fig. 1B we see the clear depletion at low velocities for powers nw1, which we aimed to achieve. The extrapolated n~1 case however does not show a significant depletion. Note that extrapolation of the fit for large powers does not leave the 0vf ƒ1 range till m §129,082. We conjecture that the marginal velocity distribution Eq.

Fits to myelinated fibre diameters in human corpus callosum
How well does the dispersive propagator and its distancedependent f jk v jr ð Þ and marginal f jk (v) velocity distributions, as well as the difference propagator and distributions derived from it, reflect physiological reality? This is a difficult question to answer since there are surprisingly few studies that have attempted to experimentally quantify the distribution of cortico-cortical conduction velocities in animals or humans. Existing experimental estimates can be divided into two groups: those based directly on conduction latencies, for which the distance-dependent velocity distribution f jk v j r ð Þ is appropriate, and those based on the transformation of histologically determined axon diameters, to which the marginal velocity distribution f jk (v) applies. Estimates of cortico-cortical conduction velocities obtained using these approaches cover a wide range, and depend on whether the fibres are myelinated or unmyelinated. For example, myelinated fibres of the corpus callosum are found to have an order of magnitude variation in diameters (0:25{2:25 mm in rat, rabbit, cat and monkey [42][43][44][45]), with conduction velocities expected to vary roughly linearly with these different calibres. Furthermore, strong regional differences can occur, for example in monkey callosal latency measurements yield a median of 7:0 ms {1 [46], whereas in visual cortex one obtains only *3:5 ms {1 [47]. In the following we will concentrate on fibre diameters and hence the marginal velocity distribution f jk (v), since here some fairly detailed data sets are available. Furthermore, the analysis of latency measurements requires knowledge about the distance between brain areas and adds uncertainties concerning the precise time when transmitted impulses actually lead to a measurable response. However, we will indicate below where latency measurements may solve ambiguities in our fits to data.
For myelinated axonal fibres conduction velocity is found to be linearly related to fibre diameter v~kd. The constant of proportionality k however is not well determined. Below we first concentrate on the work of Aboitiz et al. [31], since they provide empirical data for the distribution of callosal axonal diameters in human brains. That paper uses k~8:7 ms {1 mm {1 . But for example data summarised in Boyd and Kalu [48] suggest that for myelinated axonal fibres with diameter v10 mm the linear scale factor should be rather k~4:5{6:0 ms {1 mm {1 . However, we will see below that this uncertainty does not influence our data fit directly, but merely scales its result. Aboitiz et al. [31] obtained the number of fibres over a given threshold diameter in the corpus callosum of twenty human brains (10 males and 10 females). To this purpose saggitally sectioned and stained post-mortem callosal pieces were examined using light microscopy. In addition electron micrographs were used for one brain. A summary of their data suitable for our purposes is given in Tab. 2. Note that in this table only the last four rows and first two columns contain their directly measured data. The first row and the last two columns are estimates based on other approximate measurements also reported in [31]: On one hand we have subtracted the number of unmyelinated fibres, and on the other hand we have estimated the full unthresholded count. For details see the caption of Tab. 2.
Aboitiz et al. [31] counted the number of fibres over a given observed diameter threshold d obs . The observed diameters must be corrected for an estimated 65% tissue shrinkage due to formalin fixation and paraffin embedding [31]: d~d obs =j with j^0:65. This general shrinkage fortunately maintains the linear relation to conduction velocity: v~k=j : d obs . In order to fit this thresholded data, we calculate where N represents the number of all (myelinated) fibres in their corpus callosum sample and f jk is the marginal velocity distribution Eq. (32) of our newly proposed dispersive propagator. N thr is then the predicted number of fibres having conduction velocities larger than v. Note that thanks to the linear relationship of diameter to velocity, we can directly compare this to the experimental count N thr of the number of fibres with a diameter larger than d obs . We will then fit the optimal parameters N and d jk , and can relate the latter to the characteristic velocity as v jk~k =j : d jk . This means that the substantial uncertainties for the velocity scale factor k does not directly influence our fit. If k becomes more precisely known the new v jk can be obtained simply by multiplication. For reporting velocities we will use the factor k=j~8:7=0:65 ms {1 mm {1 in the following. An effect not covered by the general shrinkage factor j is the possibility of differential shrinkage of the tissue, i.e., fibres of different diameters may have shrunk at different rates in the preparation. Little is known about such effects. Furthermore, fibres typically have a somewhat irregular ''oval with dents'' cross section in practice, leading to uncertainties in precise determinations of the diameter. Finally, both observer error in the tedious task of counting thousands of fibres and equipment limitations (in particular for small diameters) come into play. For all these reasons it is likely that the d obs of Tab. 2 should be considered to have some error. In order to take into account all these uncertainties, in particular the unknown differential shrinkage error, we repeat the data fit four times with s dse =d obs~f 0%,2%,4%,6%g.
In Fig. 4 we show the result of fitting N and d jk for powers n from one to ten. We have repeated the fit in steps of 0.1 in order to obtain smooth curves, but as discussed above only integer powers allow easy computation in terms of PDEs. Shown is the probability of obtaining a x 2 equal to or greater than the actual x 2 , assuming that the data is drawn from the model for a selected n using best-fit parameters. This we will consider as the confidence level of the model with this particular n. We use here and throughout ''generalized chi-square-fitting'', which takes into account errors in both dependent (y) and independent (x) variables at every i-th data point using s 2 i~s 2 y,i z(Ly=Lx) 2 i s 2 x,i , to compare our predicted marginal velocity distributions with the empirically observed data. More advanced approaches, for example those based on Bayesian inference, could in principle give statistically more robust and informative estimates of model parameters. However, the kind of data available to us, from a purely practical point of view, limits the advantages one could obtain with more involved statistical analyses. On one hand, we  Light microscopy counts (first two columns) are from Tab. III in [31]. The counts for w1{5 mm used Loyez stains of only myelinated fibres, but w0:4 mm represents Holmes stains, which include unmyelinated fibres. Electron microscopy revealed ''about 16%'' unmyelinated fibres in the (three segments of the) genu and ''usually less than 5%'' in the other parts [31]. Using Fig. 1  use here aggregate data from publications, not individual observations, i.e., counts per area and per individual. On the other hand, human data is too scarce and we will see below that the rat data shows systematic deviations from the models. However, our results are sufficiently clear to exclude the longwavelength velocity distribution for all data, and motivate the use of the dispersive propagator for human and the difference propagator for rat data. In Tab. 3 we collect the results for maximum confidence level, i.e., minimum x 2 . We see that depending on whether the diameter uncertainty is larger or smaller, integer powers n~3 and n~4 are favoured, respectively. The fitted number of all myelinated fibres N remains well inside one standard deviation of the estimate N thr w0:0mm in Tab. 2, but is systematically larger and grows for larger diameter uncertainties. The fitted diameters d jk and characteristic velocities v jk are lower for larger assumed diameter uncertainties. But this mainly reflects the lower fitted powers n, since the extracted mean diameters and velocities remain similar, i.e., larger n imply ''slower'' distributions Not surprisingly, larger assumed errors allow better fits, but fit quality is generally satisfactory. Re-fitting with integer n where the best fit is obtained with non-integer n yields similar fit quality with somewhat changed parameters. Given our lack of knowledge concerning the precise diameter uncertainty, it is probably best to consider the 0% case with n~4, N~1:889 : 10 8 , v jk~1 8:74 ms {1 and the 6% case with n~3, N~1:935 : 10 8 , v jk~1 4:91 ms {1 as reasonable limiting cases. They have confidence levels of 51.41% and 68.17%, respectively. The quality of these fits is apparent in Fig. 5. Note that the 0% case predicts d mode~0 :4667 mm (v mode~6 :246 ms {1 ) and the 6% case d mode~0 :4211 mm (v mode~5 :636 ms {1 ). The difference of diameters predicted from these limiting cases is hence likely too small to be detected directly from slice measurements. However, larger n mean overall ''slower'' diameter distributions. A fit to diameter data naturally reduces the impact of n on predicted diameters, but it does so by compensating with an increase of the characteristic d jk . If conduction latencies for large distances are roughly *v jk , as speculated above, then measuring the resulting larger difference between v jk~1 4:91 ms {1 and v jk~1 8:74 ms {1 may help distinguishing the n~3 and n~4 fits experimentally. This nicely demonstrates the (speculative) complementarity of diameter and latency measurements. Conduction latencies for callosal fibres in rhesus monkey gave velocity estimates of median 7:0 ms {1 [46]. This may suggest a preference for lower v jk , i.e., n~3 and/or the lower k values of Boyd and Kalu [48], if one assumes that the inter-species difference between humans and monkeys is not too drastic. For the n~3 case with  54) was varied in steps of 0.1 for four different uncertainties of the observed threshold diameters s dse =d obs~f 0%,2%,4%,6%g. The assumed relative diameter error reflects mainly differential shrinkage. As confidence level the probability that x 2 is greater than the fitted x 2 is shown. doi:10.1371/journal.pcbi.1000653.g004  (54). Fits were repeated assuming uncertainties 0%, 2%, 4%, and 6% of the observed threshold diameters s dse =d obs for n varying from 0 to 10 in steps of 0.1. The values reported here are those of the confidence level peak, i.e., the minimum x 2 , cf. Fig. 4. Where power n was not an integer at the peak, we also provide the fit with the closest integer n, shown in square brackets. For comparison, we repeated this procedure with difference Eq. (55). Difference fits are indicated by a m~1 in square brackets, the value we used throughout to minimize the fraction of small diameter fibres, and n 1 <n, N<N tot , d 1 <d jk , v 1 <v jk for tabulation. We show results for minimal and maximal diameter uncertainties, again also constrained to integer values.
where n 2~n1 zm with n 1 ,m[N 1 . Furthermore, d 2~ffi ffiffiffi ffi n 2 p = ffiffiffiffi ffi n 1 p fd 1 and w 2~z w 1 with z and f given by Eqns. (44) and (52), respectively. We choose m~1 for the difference fits, because for this m difference and dispersive distributions are most dissimilar, whereas for m?? they become the same. We have checked that larger m indeed produce results closer to the dispersive fits. Nevertheless, even for m~1 we find confidence levels basically identical with dispersive fits of the same order, see Tab. 3. The difference fits are also shown in Fig. 5, and the similarity to the dispersive curves is evident. The only marginal improvement is that the fitted N tot are slightly closer to the experimental value for N thr w0:0mm of Tab. 2. Thus our current data for humans is too scarce and imprecise to warrant the introduction of the more complicated difference model, which requires twice the computational effort. A useful fit of the data in Tab. 2 using the long-wavelength Eq. (36) cannot be obtained, since its shape is too much at odds with what is required by the data. It is hence also not clear how to best compare results for the dispersive and the long-wavelength propagator, respectively, since their velocity distributions are so different. One possible suggestion is to match their median velocities, in which case for n~3 one obtainsṽ v jk^0 :5887 : v jk^8 :782 ms {1 , while keeping N N~N. Another possibility is to simply ignore the two data points at largest threshold diameters, which constrain the overall shape, and fit only the first three. Then a fit for s dse =d obs~6 % with a probability of 85.43% for x 2 w0:0337 can be obtained. Parameters areÑ N~1:680 : 10 8 andd d jk~1 :026 mm?ṽ v jk~1 3:73 ms {1 . Both of these possibilities are also shown in Fig. 5. We can see that the the long-wavelength propagator matches at low threshold diameters either the new propagator or the first three data points, according to our choice. For comparable numerical simulations one should also adjust the connectivity decay length according to Eq. (25), e.g., s s jk^4 :930s jk for n~3. While we strongly recommend using our new propagators, the long-wavelength one perhaps remains attractive for its computational simplicity. But in future one should then use such appropriately ''matched parameter values'' forṽ v jk and s s jk . It is interesting that the exponential propagation decay time of Eq. (6) and Eq. (11) with r?v jk t turns out to differ substantially: s s=ṽ v jk^8 :370s jk =v jk for n~3. This suggests that the dispersive propagator acts more locally than the long-wavelength propagator.

Fits to fibre diameters in rat subcortical white matter
We now turn to animals, where more comprehensive data is available. Partadiredja et al. [36] have recently provided extensive data on axon diameter distributions in rat. As mentioned above, there appears to be a general trend in lower mammals that less of the small diameter fibres are myelinated. Rat data hence provides a convenient test for our difference propagator constructed to deal with such depletion, since we would expect more small diameter fibres and hence easier fits for other animals closer to humans, for example macaques. Furthermore, unlike the human data used above and other data sets from animals, Ref. [36] resolves diameters very finely and hence allows us to pinpoint the strengths and weaknesses of our ansatz. Partadiredja et al. [36] provide their fibre count data in terms of the total densities r exp of axons per 100 mm 2 and corresponding percentage histograms p exp dependent on fibre diameters, see their Tabs. 1 and 2 and Figs. 4, 5 and 6. They provide electromicroscopic results averaged over six adult male Wistar rats, but differentiated according to myelinated and unmyelinated axons of frontal, parietal, and occipital subcortical white matter from both left and right hemispheres.
Partadiredja et al. [36] found differences between left and right hemispheres only for parietal unmyelinated fibres with appreciable statistical significance (probability ƒ4:5%). Furthermore, an independent check with a second set of data, albeit at lower magnification, did not confirm even this difference. Thus it is reasonable to average their data for left and right hemispheres. However, it remains difficult to estimate appropriately the errors on their p exp bins by only comparing data from left and right hemispheres. Considering measured mean calibres, they found only one significant regional difference (probability ƒ1:3%, parietal vs. occipital) for unmyelinated axons and one marginally significant one (probability ƒ5:6%, frontal vs. occipital) for myelinated axons. Hence we will proceed here by averaging binwise over the six p exp histograms (left and right for frontal, parietal and occipital) available each for myelinated and unmyelinated axons, and simply use the corresponding unbiased estimator of the standard deviation in our fits. It is possible that regional differences could be described with a more sophisticated procedure, but this is sufficient for a parsimonious theoretical description and judging  Tab. 3. For the dispersive propagator the n~3 and n~4 fits are shown, which are optimal assuming s dse =d obs equal to 6% and 0%, respectively. This relative diameter error (magenta error bars: 6%) reflects mainly differential shrinkage. Corresponding difference propagator fits are also shown, which have basically the same confidence levels. Thus these data cannot distinguish the dispersive and difference models, and the former is preferred for its computational simplicity. For the long-wavelength propagator a reasonable fit with Eq. (56) to all data cannot be obtained. Two curves are shown: one matching the median velocity of the dispersive n~3 case, the other fitting only the first three data points with s dse =d obs~6 %. doi:10.1371/journal.pcbi.1000653.g005 the suitability of our ansatz. We will always take direct averages of the percentage histograms instead of using weighted sums. Considering for simplicity only the average over left and right hemispheres, we thus use p exp l&r~( p exp left zp exp right )=2 instead of introducing total density weights p exp l&r~( r exp left p exp left zr exp right p exp right )= (r exp left zr exp right ). This minimizes the problem of correlated errors, since the errors on r exp were extracted from the same data as the p exp histograms. To predict densities for different diameters one should multiply the averaged p exp with the likewise averaged r exp .
We relate the marginal velocity distribution Eq. (51) once more linearly v~k : d to diameters (the diameters in [36] are already corrected for shrinkage), and obtain for a distribution in bins ½d i ,d iz1 with d i~i : Dd and i[N 0 : where n 2~n1 zm with n 1 ,m[N 1 . Furthermore, d 2~ffi ffiffiffi ffi n 2 p = ffiffiffiffi ffi n 1 p fd 1 and w 2~z w 1 with z and f given by Eqns. (44) and (52), respectively. Below we wish to compare the quality of dispersive and difference fits to rat data. The equivalent formula for the dispersive propagator from Eq. (32) is It is also easy to obtain the corresponding result for the longwavelength approximation from Eq. (36)p However, we will see that sincep p i :0 for d i wd d jk , no reasonable fit can be obtained. As in the human case, the complete lack of a large diameter (high velocity) tail in the long-wavelength distribution is at odds with experimental data from rat. The following discussion of the p tot i and P tot in Eq. (57) applies likewise to the equivalent quantities in Eqns. (58) and (59): If the fitted probability norm P tot deviates from 100%, then this indicates that the theory prefers a different total density of axons than the experimental mean, namely r tot~Ptot r exp . For comparisons with the experiment we further re-norm Eq. (57) by multiplying the predicted p tot i with a factor P tot =( P p exp i w0 p tot i ). Then the sum of the predicted p tot i over only those bins where the experimental data p exp i w0 already yields P tot . This adjusts for the systematic mismatch between the experimental data, which assigns 100% to the total as sum over those bins which have non-zero empirical entries, and the model, which assigns 100% to the total as sum over all predicted bin counts. Thus we can now truly expect P tot~1 00% from the fit.
While the linear relation v*d is widely accepted for myelinated axons [49], it is currently not clear how conduction velocity is related to diameter in unmyelinated axons. Theoretical results [50][51][52] tend to favour a v* ffiffiffi d p dependence. Experimentally one has found varying results, from squid v*d 0:61 over crab v*d 0:75 to mammalian C fibres v*d, see for example [49] and references therein. Studies of sensory neurons in cat have also suggested a linear relationship [53]. Since currently the situation  [36]. s dse =d obs~3 % (magenta error bars) has been assumed to reflect mostly differential shrinkage, but fit dependence on this is mild. Fit results using the difference Eq. (57), and its dispersive counterpart Eq. (58), are collected in Tab. 4. For unmyelinated axons the optimal fit with n 1~2 , m~1 is shown. For comparison, the optimal n~3 fit with the dispersive propagator is also displayed. It is viable, but has a three times larger x 2 . For the long-wavelength propagator a reasonable fit with Eq. (59) to all data cannot be obtained. Two curves are shown for illustration: one matching the median velocity of the difference n 1~2 , m~1 case, the other fitting only the first four data points. doi:10.1371/journal.pcbi.1000653.g006 ð57Þ is inconclusive, we use a linear relationship also for unmyelinated axons, but naturally with a lower k than for myelinated ones. As mentioned, this allows us to fit diameters directly and scale the result to velocities. If we assumed for example v~e ffiffiffi d p instead, then we would have to relate distributions nonlinearly f (v)~2vf (d~v 2 =e 2 )=e 2 . This would inconveniently turn diameter bins of the same size into different size velocity bins. Nevertheless, it is important to note that sublinear diameter powers sharpen up the velocity distribution as compared to the linear case. This is in general detrimental for the quality of our fits. Thus the linear fits to unmyelinated fibre data we provide below need be considered as the 'best case scenario'. In order to compare better with the previous fits to human data, we again use k myel:~8 :7 ms {1 mm {1 and follow Tab. IV in Aboitiz et al. [31] as well in setting k unmyel:~3 :2 ms {1 mm {1 , which is based on callosal rabbit data in Ref. [54]. We note once more that for the linear case different assumed k simply re-scale our fit results given below.
Fits of our difference propagator Eq. (57) to the data by Partadiredja et al. [36] are collected in Tab. 4, and compared there with corresponding dispersive fits using Eq. (58). The results for unmyelinated fibres are displayed in Fig. 6. The fit to unmyelinated fibre diameters has a proper optimum concerning the difference propagator model, i.e., upon trying n 1~1 , . . . ,16 and m~1, . . . ,15 we find an optimal fit for n 1~2 and m 1~1 (and thus n 2~3 ). We see that this fit is excellent with a confidence level of 99.997%, likely indicating an overestimate of the errors. Keep in mind though that this is the 'best case scenario' linear fit, with lower powers in the relation between velocity and diameter fit quality would deteriorate. Furthermore, it is noteworthy that P tot is close to 100% and that the mean diameter of the distribution corresponds very closely to the one estimated from experimental data. This further confirms that the fit performs well. However, and perhaps not surprisingly, the unmyelinated diameters can also be fit with the dispersive propagator, as shown in Tab. 4 and Fig. 6. Once more we find a proper optimum, although for n~3. All criteria for a very good fit remain: the confidence level remains high at 91.74%, P tot is close to 100% and the mean diameter is close to the experimental value. However, we see that x 2 has actually gone up by a factor of three as compared to the difference propagator fit. Due to the much higher x 2 , the quality of the dispersive fit is considerably more sensitive to the uncertainty in the relation of velocity to diameter. Inspection of the fit curves in Fig. 6 also suggests that the dispersive fit has a trend of being too wide.
The large diameter tail in the data, which precludes any direct fit with Eq. (59), is obvious in Fig. 6. One can try once more to match the median velocities of the long-wavelength approximation to that of a more viable fit. From the n 1~2 , m~1 difference fit one obtainsṽ v jk^0 :8271 : v 1^0 :6609 ms {1 . For a fair comparison, we adjust the long-wavelength data normP P optimally for the number of bins wherep p i =0. Sinced d jk^0 :2065 mm, this includes the first five bins and yieldsP P~73:95%. Alternatively one can ignore again the data points at large threshold diameters. The best fit is possible for the inclusion of the first four bins, wherẽ P P~55:56% andd d jk~0 :1656 mm (ṽ v jk~0 :5299 ms {1 ), for s dse =d obs~3 % with a confidence level of 72.52% for x 2 w0:6427. The confidence level for including the first three bins would be 51.26%, and for the first five bins 15.64%. Both the long-wavelength prediction matched in median velocity and the one fit to the first four bins are displayed in Fig. 6. Note that the long-wavelengthp p i basically rises monotonically with diameter and then suddenly drops to zero. The only slight complication arises for the last non-zerop p i , which can rise or fall as compared to the previous bin at smaller diameters. Yet an extended large diameter tail as seen in the data is impossible to achieve.  (58). Dispersive fits are indicated by a missing m value, and n<n 1 , P<P tot , d jk <d 1 , v jk <v 1 for tabulation. s dse =d obs~3 % was used for all fits, but dependence on this was mild. For myelinated fibres the difference fit had no optimal n 1 , hence several orders were tried as indicated by n 1 in square brackets. However, m~1 was optimal for any chosen n 1 . The same holds true for the dispersive fit, and matching n were tried. The entries marked with a Ã show fits made to diameters §0:2 mm only, i.e., without the first four (myelinated) data bins. Then optimal fit orders exist as shown. Here two sets of x 2 and confidence level are given: in brackets for the large diameters, without brackets compared to the full data. Unmyelinated data directly leads to the shown fits with optimal fit orders. SdT and s d are compatible with the corresponding mean over values in Tabs Turning now to our fit results for myelinated fibres, see Tab. 4, we find that the difference propagator has some trouble matching the data. In Fig. 7 this is illustrated by two curves for n 1~4 and n 1~8 , respectively, with m~1 in both cases. Essentially, the experimental distribution is more sharply peaked around 0:3 mm than the difference propagator can easily accommodate. Since the difference propagator becomes more sharply peaked for larger n 1 , higher powers always provide a better fit in the tested range n 1 ƒ16, i.e., we cannot find an optimal n 1 for the model. However, for any given n 1 one finds that m~1 is optimal, since that maximizes the skewness of the distribution. That P tot is only about 75% also indicates that our fit has trouble matching the sharp maximum. That said, formally one finds reasonable confidence levels for higher powers of n 1 , e.g., 53.45% for n 1~8 . Furthermore, the mean diameters of the distributions are well compatible with the experimental value. While acknowledging the difficulties, we hence conclude that the difference propagator is sufficient for a rough fit even to rat data. Fitting the longwavelength propagator to these data is of course hopeless, due to the extended large diameter tail. Since we have already considered artificial matching procedures in the unmyelinated case, we do not discuss any long-wavelength fits here. The dispersive propagator also prefers large n without proper optimum. But if we fix n to the same value as the n 1 of the difference propagator, then we find roughly a two times larger x 2 in the dispersive case. This then implies negligible confidence levels for the fit, i.e., the dispersive propagator can be considered as excluded for the myelinated rat data. We show in Fig. 7 two corresponding dispersive curves with n~4 and n~8, respectively. It is obvious that compared to their difference counterparts they are primarily less able to accommodate the sharp peak around 0:3 mm.
These results can be summarized also as follows: the depletion of small diameters fibres in the experimental data appears to be even stronger than predicted by the difference model, and excludes the dispersive model. In order to demonstrate that the small diameter data is the culprit, we have repeated the fits, but removed the small diameter bins one by one. We find that after removing the first four bins, and thus for considering only diameters greater than 0:2 mm, both the difference and the dispersive fit acquire optimal fit orders, namely n 1~6 , m~1 and n~7, respectively. These fits for larger diameters are also shown in Fig. 7, and are indicated by a Ã in the legend. As one can see in the figure and in Tab. 4, fit quality is excellent for large diameters for both models. But if one uses the so obtained parameters and compares to the full data set including the small diameter bins, then the confidence levels become negligible. Though again the x 2 of the dispersive model is about two times larger. It is possible that some experimental problem exists that leads to a systematic underestimate of the number of small diameter myelinated fibres, though we are not in fact aware of any. If that were the case, then the large diameter fits might be closer to reality. Furthermore, the large diameter fit for the difference propagator is actually in accord with the two smallest diameter bins. Hence one could use it instead of say the regular n 1~8 fit, in order to trade a mismatch in the third and fourth bin for an improved description of the peak. Myelinated diameter counts for higher mammals, which show less depletion at small diameters, should be described more easily with the difference propagator. Indeed, this is also suggested by the Figure 7. Fits to binned counts of myelinated fibre diameters in rat subcortical white matter. Data and fits are obtained as for Fig. 6, but using the myelinated counts. Two regular difference fit curves are shown: n 1~8 and n 1~4 , with m~1 in both cases. Systematic deviations from data around 0:3 mm are obvious, but fit quality remains tolerable with a confidence level of 53.45% for n 1~8 . Even larger n 1 can increase the confidence level to about 70%. For comparison, dispersive fits with orders n~8 and n~4 are also shown. Their x 2 is almost a factor two larger, rendering their confidence level negligible. Fits with the long-wavelength propagator are not show, but fail drastically, cf. Fig. 6. The curves marked with a Ã show additional fits for diameters §0:2 mm only, i.e., without the first four data bins. Then one can find optimal fit orders for both propagators. These fits are of comparable, excellent quality compared to the reduced data set. But both predict too many small diameter fibres, and hence have negligible confidence levels compared to the full data set, with the dispersive x 2 again being about two times larger. doi:10.1371/journal.pcbi.1000653.g007 success of our own fit with the dispersive, and hence ''nondepleted'', propagator to human data.
We can now use the fit to unmyelinated rat data to speculate about the human case, for which we have not enough data available for an independent fit. Let us assume that like unmyelinated rat subcortical white matter, also human unmyelinated callosal fibres can be fit with a n~3 dispersive propagator. For human callosal fibres, we see from Tab. 2 that 9.19% of fibres are unmyelinated, whereas for rat subcortical white matter we derive from the mean numbers in Tabs. 1 and 2 of [36] that 83.35% of fibres are unmyelinated. We can use these fractions to construct combined marginal velocity distributions in order to understand overall activity conduction properties: where we have used the labels of Tab. 5 as subscripts to indicate alternatives. In Fig. 8 we show these combined distributions, and the respective myelinated and unmyelinated contributions. To disentangle the curves vf (v) rather than f (v) is shown. Thus the area of these curves is normed to mean velocities, rather than to one. Furthermore, to give some feeling for the remaining uncertainty even in our ''best fits'', we show bands using the minimum and maximum envelopes of Eqns. (61) and (62). Thus for example the lower border of the ''rat -myelinated'' band is computed as 0:1665 : There are of course considerable caveats: the human unmyelinated part is derived speculatively, the rat myelinated part is only a rough fit, and the unmyelinated estimates are in general plagued by the uncertain relation of diameter to velocity. Nevertheless, we expect that the clear differences one can observe here will hold true at least qualitatively: In rat subcortical white matter there are two modes, a dominant, sharp one at low velocities and a broad one at higher velocities. In human corpus callosum one finds only a single, very broad mode at high velocities. One would have to lower the ratio of the myelinated to unmyelinated k from 8.7/3.2 = 2.7 to about 1.9 to turn the second rat mode into a high velocity shoulder, and further down to about 1.0 to obtain a smooth unimodal distribution. It is biologically implausible to assume that the k ratio could be so low, since that would abandon the distinction between fibre types. Furthermore, a sublinear relation of velocity to diameter in rat would sharpen the distinction between the modes even more. It is hence likely that rat subcortical white matter operates in two distinguishable velocity regimes, whereas human corpus callosum features only a single one.

Turing instability analysis
Following Coombes et al. [30] we investigate the consequences of the new dispersive propagator in terms of a Turing instability analysis. The Turing instability analysis represents the standard approach to understanding the emergence of spatio-temporal Table 5. Summary of best propagator fits recommended for use with human and rat data.  patterns of activity in spatially continuous non-linear dynamical systems [55][56][57][58]. Specifically it enables the determination of the conditions under which the stability of a homogeneous steady state is lost and the types of patterns of activity that subsequently emerge. This analysis method has been of great utility in understanding self organized pattern formation in a range of physical, chemical and biological systems. In order to facilitate comparison with previously developed long range propagators, we explore the stability of the homogeneous steady state for a Wilson-Cowan or Amari style neural field model in which the mean soma membrane potential is given by Here g jk (t) corresponds to the time course of a unitary postsynaptic potential (PSP) and w jk (x,t) represents the total rate of arrival of presynaptic impulses to neuronal population k arising from neural population j. We choose the bi-exponential to model PSPs. Therefore the system of equations for the Turing instability analysis are where k k is a gain parameter in the firing rate sigmoid S k . The homogeneous steady state h k (x,t)~h Ã k is then given by h Ã k~Pj w 0 jk S j ½h Ã j zh 0 k . A linearization around this state with h k (x,t)~h Ã k zh k e lt e ik : x and like perturbations of u jk and w jk yields the system of equations where S' j ½h~LS j =Lh and c j :S' j ½h Ã j . See Refs. [55][56][57][58] for further detail on this linearization method and its application to the Turing instability analysis.
Nontrivial solutions for h k will only exist for det½D(k,l){1 :E(k,l)~0, where 1 is the identity matrix of appropriate dimension. Solutions to E(k,l)~0 then yield a continuous spectrum of eigenvalues, l~l(k), that define the dispersion relationship. Clearly each spatial mode k will be stable if the real parts of the corresponding eigenvalues l(k)v0, and thus the homogeneous state will be stable to all perturbations if <l(k)v0 for all k. As various model parameters are changed, we expect a critical point will be reached for some k~k c where the real part of the corresponding eigenvalue l(k c ) becomes zero. By parametrically moving beyond this critical point the eigenmodes having critical wavenumber k c and critical frequency v c :=l(k c ) will start to grow, leading to the emergence of spatio-temporal patterns of activity. The expected type of emergent activity can be inferred from the values of k c and v c . If k c~0 and v c =0 then we expect to see the emergence of spatially uniform periodic oscillations. If k c =0 then we expect to see the emergence of spatial patterns of activity that can either be periodic in space but constant in time (v c~0 ), or periodic in space and time (v c =0). These three bifurcation scenarios are typically referred to as Hopf (k c~0 , v c =0), Turing (k c =0, v c~0 ), and Turing-Hopf (k c =0, v c =0) bifurcations, respectively.
For computational purposes it is preferable to split E(k,l)~0 into real and imaginary parts and define l:nziv: where q~fa jk ,b jk ,n,w 0 jk ,v jk ,s jk ,k k g indicates the chosen set of model parameters. Solutions to Eqns. (70) and (71) for a given set of parameters q thus yield curves n(k,v; q) in the (n,v)-plane parameterised by k, see for example the insets in Fig. 9B. Formally a Hopf bifurcation occurs when k c~0 and v c =0 which from Eqns. (70) and (71) gives the condition A Turing-Hopf bifurcation occurs when k c =0 and v c =0 and requires that the solution trajectory n(k,v; q) should be a tangent to n~0 at k c ,v c ð Þ, i.e., in addition to k c and v c satisfying Eq. (73) can be derived by noting that the total derivatives of Eqns.
(70) and (71) with respect to k are by the chain rule Solving these two equations by eliminating dv=dk yields as conditions for the bifurcation. The derivation of Eq. (77) proceeds in a similar manner to that of Eq. (73). In principle one also has to check that such tangential solutions are locally rightbounded, with the local turning point being least stable, but for the dispersive, difference, and long-wavelength propagators under consideration we have found this to be always the case in practice.
In the following, we will investigate the existence of these bifurcations by changing one model parameter q a and solving the equations for v c , k c , and a second model parameter q b . As in Coombes et al. [30], for simplicity we consider only two neuronal populations: excitatory (k~e) and inhibitory (k~i). We further simplify the PSP time courses by setting both a ik~bik~1 and a ek~bek~1 . Explicit synaptic delays are not modelled. Because in neocortex excitatory connections have much greater lateral extent than inhibitory connections, we assume that s ek ws ik , and here set s ik~1 and s ek~2 . Connectivity weights represent local dominance of inhibition with w 0 ek~1 and w 0 ik~{ 4, respectively. Uniform axonal conduction velocities v jk~v and firing rate functions S j~S are assumed for simplicity. For subsequent numerical simulations, and without loss of generality, we set h Ã k~0 so that the linearized gain c~k=4. Fig. 9 shows the results of the Turing instability analysis for the dispersive and long wavelength propagator models. Figure 9A shows the critical curves in the (v,c) plane. Above each of the respective curves a homogeneous steady state succumbs to dynamical instabilities for k c~0 (bulk oscillations, Hopf) and k c =0 (travelling waves, Turing-Hopf), with the lower curve determining the actual bifurcation at a given v. Neither model gives rise to a Turing bifurcation within the admissible parameter space, i.e., Turing bifurcations only occur for negative c. This is in contrast to Steyn- Ross et al. [59], in which the effects of both chemical and electrical (gap junction) synaptic transmission are modelled in a mean field model that includes a long-wavelength propagator. They observed stationary Turing instabilities when homogeneous driving terms similar to h 0 k were varied. Both the dispersive and long wavelength propagators are capable of exhibiting Hopf and Turing-Hopf instabilities. As Fig. 9B shows for the chosen parameter sets, the Turing-Hopf critical wavenumber k c exists only above a finite velocity v for both dispersive and long-wavelength models. In all cases the k c rises smoothly with v and asymptotes to a maximum value. Considering the critical wavelength l c~2 p=k c as a typical size, increasing v will hence change bulk oscillations first into large travelling waves, which then contract to some minimal size. While the boundaries of instability are of similar shape in the long-wavelength and dispersive propagators, there are nevertheless differences between the loci of the two curves that may be of considerable physiological relevance. As can be seen in Fig. 9A, loss of stability occurs for significantly smaller values of the linearized gain c in the dispersive model with nw1 as compared to the long-wavelength one, for a given characteristic conduction velocity v. In general this remains true for a reparametrization of the bifurcation curves in terms of the mean, median, or mode of the corresponding velocity distributions, cf. Tab. 1. Because the linearized gain c:S'½h Ã ~kS Ã (1{S Ã ) where S Ã~S ½h Ã is the steady state firing rate, cf. Eq. (66), changes in c can be achieved by alterations in the steady state firing rate via Dc~kDS Ã (1{2S Ã base ). Thus for the nw1 dispersive propagators a smaller change DS Ã from a given basal firing rate S Ã base is required to induce pattern formation, as compared to the long-wavelength propagator. This may follow more closely the biological situation, where a range of metabolic and energetic constraints need to be negotiated.
In Fig. 9C we show the critical frequency v c only of the less stable bifurcation, which actually determines the instability. The v c of the propagators are seen to transit smoothly from Hopf to Turing-Hopf with increasing velocity. However, the long-wavelength v c increases more quickly with velocity than the dispersive ones, except for the n~1 case which is a close match. Thus at a given velocity, nw1 dispersive travelling patterns will emerge at lower critical frequencies than long-wavelength ones. Fig. 9D displays the critical phase velocity c~v c =k c of the emerging patterns. We find a lower v limit for which c formally diverges, since k c ?0 in this limit. Both the dispersive and the longwavelength critical phase velocities then rise mildly for larger velocities. It is known that developmental changes to the diameters and myelination of axonal fibres occur and partly depend on activity feedback, see for example [60] and references therein. Although highly speculative, it is conceivable that the transitions between bulk oscillations and travelling waves in response to changing conduction velocities, which we have just discussed, could provide a relevant feedback mechanism. That the phase velocity c remains close to independent of v above a thresholdparticularly so for larger n dispersive propagators, less so for smaller n and the long-wavelength case -may then be significant for connectivity development. Obviously significant differences exist between the dispersive and long-wavelength models, especially for larger n. Such differences likely also occur in the bifurcation structure of other parameter planes. The biological implications of these dynamical differences require future detailed investigations, which need to go beyond our qualitative considerations here by restricting the parameters more specifically to experimentally allowed ranges.
In order to test the predictions of our linear stability analysis we have performed numerical simulations of Eqns. (66) to (68) over suitably chosen domains. Tab. 6 shows the results of comparisons between the spatio-temporal properties of the numerical simulations for (v,c) just beyond the Turing-Hopf bifurcation, and the corresponding linear predictions. As can be seen there is excellent agreement for a range of parameters and dispersive propagator orders. In all cases parallel moving stripes were seen beyond the Turing-Hopf bifurcation when integrations were continued for a long enough time (results not shown), but a range of other patterns also occurred depending on the initial conditions. Thus this system likely possesses multiple attractors. By moving further away from the Turing-Hopf bifurcation boundary more complicated, and arguably biologically more plausible, self-organizing behaviour is seen. One such example is shown in Figure 10, see also the corresponding supplementary animation S1. No attempt was made to determine whether the Turing-Hopf bifurcations were subcritical or supercritical in character, though in principle this could have been established by brute force numerical simulation or more elegantly using the method of harmonic balance [61].
Performing a Turing instability analysis for either of the difference propagators of Eq. (62) in the (v 1 ,c) parameter plane revealed a qualitative match with the corresponding n~n 1 dispersive propagator. In particular, root-loci parameterised with respect to k, see insets of Fig. 9B for representative examples, reveal that the effect of depleting low velocity fibres in accord with the rat data is to alter the most weakly damped branch of the dispersion relationship only for large values of k, with the low critical wavenumbers observed for the dispersive propagator remaining essentially unchanged. Since the critical curves for the difference propagator would basically reproduce those of the dispersive propagator in Fig. 9, we do not show our additional results for the difference propagator. These predictions also have been verified by numerically integrating Eqns. (66) and (67) using either the un-myelinated or myelinated difference propagator of Eq. (62). Thus on the basis of a linear instability analysis in the context of our current, highly simplified, neural field theory there  Fig. 9. For numerical simulations, a c sim somewhat larger than c lin c was chosen. The space-averaged 1D temporal Fourier spectrum u ee (v) was used to estimate v sim c as the maximum of ju ee (v)j. The time-averaged 2D spatial Fourier transform u ee (k x ,k y ) was used to obtain two estimates: k sim c1 as the k~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k 2 x zk 2 y q for which ju ee (k x ,k y )j is maximal; and k sim c2 as the k for which the mean of ju ee (k x ,k y )j over a circle around the origin with radius k is maximal. For the estimates grid time series of 50 time units with dt~0:01 (500 samples total) were recorded, after initial ''transients'' of 100 ( { 3000) time units were discarded. The spatial grid was 128|128 ( Ã 180|180) with discretization steps Dx~Dy~1. doi:10.1371/journal.pcbi.1000653.t006 appear to be no essential dynamical differences between the dispersive and difference propagators concerning bifurcations. However, the dynamical consequences of these propagator forms need to be investigated in future with physiologically more realistic theories of brain electrorhythmogenesis, such as those of Liley et al. [26] and Robinson et al. [25]. Furthermore, the combination of myelinated and unmyelinated fibre systems, see Eqns. (62) and (61), as well as Fig. 8, reveals differences that are likely to be significant dynamically. However, in order to make meaningful inter-species comparisons of brain dynamics, on one hand more experimental connectivity data, in particular human, is required and on the other hand local brain activity descriptions by mean field theories will have to be adjusted for different species as well.

Discussion
Understanding the physiological basis of brain dynamics requires one to account for the activity of distributed populations of cortical neurons. Modelling the details of their ongoing communication will generally form an important component of any theoretical description. Continuum mean field models (MFMs) of neural population activity [1,2,4,23,25,26] are a particularly useful theoretical tool for bridging the gap between the macro-to mesoscopic assays associated with non-invasive neuroimaging (e.g., EEG or fMRI BOLD) and our knowledge of the underlying microscopic anatomy, physiology and pharmacology. However, MFMs have faced one particularly significant technical challenge: the biologically plausible, yet computationally tractable, propagation of neuronal activity via long-range (corticocortical) connectivity. Most current MFMs have followed the pioneering work by Jirsa and Haken [27], which relied on a number of simplifying assumptions in order to derive a numerically efficient and analytically tractable ''long-wavelength'' propagation PDE. However, in doing so a substantial degree of biological fidelity has been lost, the most crucial of which involves the distribution of axonal conduction velocities. As we have demonstrated here, these PDE formulations have assumed a sharply peaked velocity distribution with a definite cut-off, a feature which is completely at odds with the available empirical evidence that instead suggests rather broad distributions. On this backdrop we have introduced two new long range propagators, the dispersive propagator and the difference propagator derived from it, which retain all the advantages of a PDE formulation but produce broad velocity distributions in keeping with the experimental measurements.
We have provided an extensive analysis of the mathematical properties of these new propagators, and contrasted them with the commonly used long-wavelength model. Of particular note are the following results: First, we can distinguish between the distancedependent and the marginal velocity distribution. The former is appropriate for the description of experimental measurements of conduction latencies, the latter can be related to the histological determination of fibre diameters in slices. Second, our new propagators predict that more distant brain areas are generally connected by faster fibres. This could be relevant for isochronicity in the brain, see for example [60]. In contrast, the long-wavelength propagator assumes essentially one conduction velocity irrespective of fibre length. Third, if conduction velocities of fibres indeed depend on distance, then typical velocities as extracted from latency and diameter measurements, respectively, are expected to differ. We are not aware that this effect has been described or systematically studied so far. For our new propagators, we can speculate that measuring activity delays over large distances should results in faster velocity estimates than deriving them from the diameters observed in local slices. It may even become possible to falsify propagator models using the constraints from these different types of data, though it is at present unclear whether the current  Fig. 9A. Spatial derivatives were approximated using finite differences on a regular square grid of 128|128 with spacing Dx~Dy~1. The resulting system of equations was rewritten as a first-order system and integrated using ode45 in MATLAB starting from random initial conditions in u ee . See also the supplementary Video S1 for the corresponding animation. doi:10.1371/journal.pcbi.1000653.g010 modelling of brain anatomy is accurate enough to allow such conclusions. Fourth, we have shown that two dispersive propagators can be subtracted from each other such that the resulting difference propagator does not exhibit unphysiological properties. This approach, which we have introduced here to deplete the number of low velocity fibres for our rat data fits, can be generalized to the construction of other, more complicated propagators as need arises.
The empirical relevance of our proposed dispersive propagator was illustrated by fitting the associated marginal axonal velocity distribution to histological measurements of axonal fibre diameters obtained from human corpus callosum by Aboitiz et al. [31], see Fig. 5 and Tab. 3. A similar fit with the long-wavelength model was simply impossible, since its functional form entirely mismatched the data. Thus this fit provides for the first time a realistic description of activity propagation in the human brain in the context of MFMs, though unfortunately only callosal data is available in the human case. In order to obtain data from other subcortical matter we turned to the extensive data of Partadiredja et al. [36] for rat. Here we had to introduce the difference propagator to account for the depletion of low velocity axons relative to human in lower mammals. The difference marginal velocity distributions were then shown to fit reasonably well the empirically derived distributions of rat axonal conduction velocities, see Figs. 6 and Fig. 7, as well as Tab. 4. Some systematic deviations between theory and experiment were however visible, caused by an even stronger low velocity axon depletion in the data as compared to our theory. However, it is known that this depletion is the less severe the phylogenetically higher the animal. Indeed we have described the human callosal data successfully here with the ''non-depleted'' dispersive propagator. Hence it is likely that rat data is a kind of ''worst case'', and human data a kind of ''best case'', for our new propagators, and reasonable fits were obtained for both.
The results of these fits further allowed us to speculate that the overall velocity distributions of rat subcortical and human callosal fibres are qualitatively quite distinct, see Fig. 8. In rat subcortical matter one finds two modes: a narrow low velocity one corresponding to unmyelinated fibres and a broad high velocity one corresponding to myelinated fibres. Whereas in human corpus callosum there is a single very broad mode at high velocities supported by both fibre types. Therefore rat data (and possibly other mammalian data) may be misleading for the purposes of MFM parameterisation in the context of understanding the spatiotemporal dynamics of human long-range connectivity. Nevertheless, caution must be exercised regarding the actual fitted parameter values due to limitations in the available data: in all cases only data aggregated across individuals and (sub-)regions was used, and this aggregate data was scarce for human and proved difficult to fit for rat. While the advantage of incorporating a high velocity tail with the dispersive and difference distributions is compelling for all data, and the depletion of low velocity fibres with the difference one is important for data from lower mammals, more robust estimates of the fitted parameters will be essential to obtain greater biologically fidelity in future MFM studies. This will depend on the availability of more and ''purer'' empirical data, as well as the use of more advanced inferential methods for the parameter estimation. For example, an analysis could be made using the Bayesian inferential framework, whereby prior beliefs about the parameters are updated using the available data. One can then simulate from the resulting posterior distributions of the parameters using Markov Chain Monte Carlo (MCMC) methodology [62]. This would have the advantage of more fully taking into account parameter uncertainty, and would allow direct probability statements to be made about the parameters. It would however introduce additional complications in terms of the model fitting process. The outcome would also be dependent on prior beliefs; in the absence of prior beliefs, one could specify prior distributions that are uninformative, but there are several ways of doing so. The frequentist approach that we have followed here is relatively simple to apply, is objective, and comparison of fitted models is straightforward.
In order to investigate the dynamical consequences of our newly proposed propagators, we have embedded them in a Wilson-Cowan or Amari style neural field formulation of local activity which while somewhat simplistic and abstract, is nevertheless more amenable to analytic treatment than biologically more realistic models such as those of Liley et al. [26] and Robinson et al. [25]. This also facilitates comparisons with previous works [30] and highlights contributions to the observed dynamics from activity propagation. Turing instability analyses were then used to characterise how spatially homogeneous steady states lose stability, and in particular how the patterns of emergent spatio-temporal activity vary as model parameters are changed. These analytic results are based on a systematic linearization of the model, but were confirmed for several cases with numerical simulations of the full equations. The difference propagator was seen to result in essentially the same bifurcation dynamics as the dispersive propagator, at least in this setting. However, considerable differences in the bifurcation dynamics between the dispersive and the long-wavelength propagator were found. Both models predict a transition for increasing axonal conduction velocities from bulk oscillation to travelling waves as dominant instabilities. But the dispersive propagator more easily transits from a homogeneous stable state to self-sustained spatio-temporal patterns. In particular it is found that pattern formation can be induced for smaller changes in neuronal firing rates with the dispersive propagator compared to the more standard long-wave length propagator for given axonal conduction velocities. The biological implications of these features are at present unclear, though it might be speculated that this represents better the biological situation, where a range of metabolic and energetic constraints need to be negotiated in order to ensure that pattern formation, and thus perception, occurs in optimal circumstances. However, an important qualification needs to be attached to these results. The emergence of self-sustaining spatiotemporal patterns of activity was predicted and simulated with isotropic and homogeneous connectivity. We did not explore here more realistic synaptic footprints, since their inclusion would have considerably complicated the qualitative picture we wished to paint by requiring mappings to coupled PDEs and multiple cortical patches [30,37,38]. Furthermore, for quantitative predictions independent connectivity data at the level of detail appropriate for MFMs is still lacking, e.g., what fraction of synapses on a local MFM neuron are associated with input from a specific distant region of cortex is generally not known with adequate precision. However, Jirsa and Kelso [63] have elegantly shown that the stability of spatial patterns can be changed by systematically varying the underlying connection topology. Even relatively simple MFMs can then undergo a series of spatiotemporal bifurcations. Since the stability of spatial patterns could also critically depend on heterogeneous connections, considerable uncertainty remains concerning the effects of conduction velocity distributions which we have reported here. It is essential that further work is performed to systematically assess the effects of conduction velocity distributions together with that of heterogeneous connectivity. In this regard it is fortunate that recent advances in modelling the latter [30,37,38] can be straightforwardly combined with our work here by changing the underlying conduction PDEs to our dispersive or difference propagators. Future studies with biologically realistic MFMs need to consider carefully whether changing the propagation model would significantly alter their predictions.
Studying the dynamical consequences of the dispersive propagator with more realistic MFMs of brain activity for example may provide greater insight into the role variations in axonal conduction velocity have in health and disease. For instance it has been hypothesised, on the basis of physiological measurement, that general anaesthetics may alter cognitive function through their effects on axonal conduction velocity [64]. A variety of general anaesthetic agents can cause increases in axonal conduction velocity of 10-20% in the peripheral nerves of human volunteers [65]. It is therefore reasonable to speculate that similar changes will occur in other myelinated axons, such as myelinated corticocortical fibres. However, more recent studies involving hippocampal tissue slices have shown no effect of the volatile anaesthetic halothane on the conduction velocities of myelinated Schaffer collaterals [66]. Our newly developed propagator, in the context of a realistic mean field theory of electrocortical activity, may help resolve the role that changes in conduction velocity have in determining anaesthetic action. Now that our novel propagators allow reasonable fits to experimental data in animal and human, we hope for a surge in theoretical investigations of conduction effects, which in turn should stimulate more targeted experimental measurements. In particular, MFMs can now include realistic activity conduction on an empirical basis in the computationally convenient fashion of a PDE for the first time. Furthermore, we expect that our result that fibre diameters and activity latencies estimate different, complementary, aspects of conduction in the brain to be a general feature of underlying velocity distributions, and hence to be of general interest beyond the specific scope of our current work. Finally, our finding that rat subcortical and human callosal fibre systems differ significantly in their velocity distributions beyond simple scaling, while admittedly speculative and clearly limited due to the comparison of different anatomical regions, is of great significance in terms of the inferences we can make about human brain activity from animal models. In particular more attention must be paid to the possible confounding effects that models parameterised on the basis of animal data have in theoretically characterising and accounting for the propagation of axonal activity in human brains.

Supporting Information
Video S1 Spatiotemporal patterns of activity produced by a Wilson-Cowan or Amari style neural field model with the dispersive propagator. The video shows a numerical simulation (1000 frames at a resolution of 0.01 time units) of Eqs.