Towards a model-based integration of co-registered electroencephalography/functional magnetic resonance imaging data with realistic neural population meshes

Brain activity can be measured with several non-invasive neuroimaging modalities, but each modality has inherent limitations with respect to resolution, contrast and interpretability. It is hoped that multimodal integration will address these limitations by using the complementary features of already available data. However, purely statistical integration can prove problematic owing to the disparate signal sources. As an alternative, we propose here an advanced neural population model implemented on an anatomically sound cortical mesh with freely adjustable connectivity, which features proper signal expression through a realistic head model for the electroencephalogram (EEG), as well as a haemodynamic model for functional magnetic resonance imaging based on blood oxygen level dependent contrast (fMRI BOLD). It hence allows simultaneous and realistic predictions of EEG and fMRI BOLD from the same underlying model of neural activity. As proof of principle, we investigate here the influence on simulated brain activity of strengthening visual connectivity. In the future we plan to fit multimodal data with this neural population model. This promises novel, model-based insights into the brain's activity in sleep, rest and task conditions.


Introduction
Non-invasive recording of human brain activity has a long history, beginning with the electroencephalogram (EEG) [1]. The EEG remains prominent both in research and in clinical practice [2] owing to its excellent time resolution, which allows, for example, the tracking of evoked potentials. In the meantime, functional magnetic resonance imaging (fMRI) based on blood oxygen level dependent (BOLD) contrast has become a standard for researching cognition [3,4], largely because fMRI BOLD can locate brain activity with millimetre accuracy. However, progress in neuroimaging is slowing, owing to fundamental restrictions in acquisition methods. For example, volume conduction limits the spatial resolution of EEG to the centimetre range, whereas fMRI BOLD relies on vascular changes with latencies of half a second or more.
A promising way forward lies in combining already available neuroimaging modalities [5][6][7][8][9]. Each modality provides a particular and distinct representation of the brain's state via its specific signal sources. For example, EEG relies on electrical and fMRI BOLD on haemodynamic sources, which relate to different aspects of neural activity, as we shall see. Furthermore, data from different modalities often have complementary characteristics. This is the case for EEG and fMRI BOLD concerning spatiotemporal resolution, as discussed above. Finally, it is possible to record EEG and fMRI BOLD simultaneously [10]. This avoids the question whether data obtained in different sessions really refer to the same brain state. Hence EEG and fMRI BOLD present a convenient test case for developing multimodal approaches. An example for the 'added value' that simultaneous EEG/fMRI can provide is the connection of fMRI resting-state networks to EEG cortical microstates recently discovered by Britz et al. [11] and Musso et al. [12] (cf. the commentaries by Laufs [13] and Lehmann [14]).
There are three basic approaches to multimodal integration [15], which are all 'model-based' in some sense. When using converging evidence, a researcher combines data argumentatively against a backdrop of established expert opinion. Clearly, such an implicit 'human mind' model is powerful, but qualitative and idiosyncratic. Data fusion combines the various recorded data directly using statistical methods. This is quantitative and repeatable, but implies some model of relations between signal sources. Such relations may be less simple than commonly assumed (e.g. [16]). Computational modelling makes explicit prior knowledge and assumptions through the process of model specification. In principle, this allows the fully objective assessment of theory in terms of a model fit to data. In practice, realistic models are often too complex for a comprehensive validation or unequivocal falsification. Nevertheless, the explication of theory through a computational model generally allows more rigorous testing.
Our NPM has been introduced previously in Bojak et al. [42]. Here, we add considerable technical detail necessary for implementing such a model. The paper is organized as follows. The following section explains how we extract the head model from structural MRI. Section 3 explains our NPM and signal expression. Section 4 presents new results for variations of specific connectivity strength. We then conclude with a discussion and outlook.

The head model (a) Surfaces extracted from structural MRI
Surface approximations for the interfaces between grey matter (GM), cerebrospinal fluid (CSF) and white matter (WM) were obtained using the CIVET software pipeline [43]. This involves a series of processing steps. Firstly, field nonuniformity artefacts are removed from T 1 -weighted structural images with the N3 algorithm [44]. T 1 -weighting is a standard MRI acquisition protocol that provides optimal intensity contrast between the tissue types of interest here. Secondly, the corrected images are normalized to stereotaxic space, and subsequently skullstripped and classified into GM, WM and CSF. Thirdly, the GM/WM interface is constructed by deforming a spherical surface mesh subject to optimization constraints. Fourthly, Laplacian GM fields are computed, and CSF skeletons are constructed in deep sulci, as guides for an expansion of the GM/WM surface to the GM/CSF interface. The resulting meshes represent the borders of cortical GM, and we use an intermediate surface for modelling; see figure 1a.
To estimate volume conduction for EEG predictions, it is also necessary to obtain skull and scalp surfaces. Using the same images, and the obtained cortical surfaces as a starting point, T 1 intensities can be sampled along rays extending outwards until a voxel containing 'air' is found, indicating the edge of the scalp. In such a T 1 -weighted image, skull tissue produces little signal, and thus its inner and outer boundaries can be determined as the edges of a 'dip' along the ray. Figure 1b illustrates this method. More details will be provided in a forthcoming publication.

(b) Pruning the cortical mesh
The CIVET mesh consists of 81 920 triangles with 40 962 vertices per hemisphere, but cortical folding can be represented faithfully with much fewer vertices. This is essential to reduce the NPM computation time, which scales linearly with the number of vertices. We will see below that 'background connectivity' scales roughly with the surface area, hence with the square of the number of vertices, and that longer edges between vertices reduce significantly the overall data transfer during the simulation. Thus, it is particularly advantageous to prune short edges. Furthermore, inspection of the CIVET mesh shows that many small triangles are 'wasted' on relatively flat parts. In other places there appear unnatural ripples, typical evidence for 'over-fitting'. These issues also are dealt with by removing short edges. Hence in each iteration of our algorithm, the shortest edge of the mesh is found. Subsequently, it is removed in a manner that we will detail next. The two triangles that previously shared the shortest edge are removed, and their other two sides collapsed to a single edge terminating in a new vertex; see figure 2a. If the new vertex were to be positioned halfway on the former edge, the curvature of the cortex would be poorly preserved. Instead, the new vertex is positioned on a circle segment that approximates the surface (cf. figure 2b): p 1 and p 2 are the vertices to be removed, and n 1 and n 2 are the surface normals at these vertices, which are defined as the average over the area-weighted normals of the triangles they belong to. First p 0 = 1 2 (p 1 + p 2 ) and n 0 = 1 2 (n 1 + n 2 ) are computed. Next, n 1 and n 2 are projected onto the plane through p 0 that is spanned by n 0 and p 1 − p 2 , resulting in n 1 and n 2 , respectively. Point c is the intersection of the lines through p 1 with direction n 1 and through p 2 with direction n 2 . Then d k is defined as the distance between c and p k for k = 1, 2, and The method fails for cases like edge 1-2 in figure 2c, since removing it would collapse two triangles onto each other. We refer to this as the 'tetrahedron', since the configuration resembles one if vertex 4 is not in the same plane, cf. figure 2d. Note that the base triangle 1-2-3 is empty, i.e. not itself part of the cortical surface. Tetrahedra are detected and removed by replacing them with single triangles; see figure 2d. We use a linear factor 0 ≤ l < 1, where l = 0 indicates the edge base points 1, 2 or 3, and l = 1 the peak 4. Thus, the tetrahedron is replaced by the base triangle for l = 0 and with lifted ones for 0 < l < 1. While l > 0 compensates better for loss of volume, it distorts the surrounding triangles. Figure 2e shows an original CIVET mesh (left) pruned to a minimum edge length of 2.5 mm (right). We have used here b = 1 and l = 0. The pruned surface has 32 408 triangles with 17 208 vertices. The overall shape and surface area remain well preserved.

(c) Parcellation and specific connectivity
Tractography based on diffusion MRI is popular for determining connectivity [45], but has significant drawbacks. Firstly, it is biased towards short-range connections and has problems where tracts are densely packed. Secondly, it does not determine the direction of activity propagation. Thirdly, it cannot find precise termination points in the GM. Data obtained through histological tract tracing methods are free of these problems, but only available from animal studies. We use here a connectivity matrix for macaque monkey (available on request), based on tracer data in the COCOMAC database [46], together with a 'regional map' (RM). This RM is a parcellation of cortex, which is sufficiently generalized to accommodate anatomical homology across primate species and uses area names that are widely recognized and convenient to use [47]. Since the CIVET pipeline is designed to obtain optimal correspondence between individual vertices of meshes extracted from different human brains, it is sufficient to map macaque brain regions to a CIVET mesh once. The assignment of vertex number to equivalent brain region stays the same for all CIVET meshes. The RM was manually delineated onto a template macaque cortical surface (F99-UA1) and, by using two landmark-based deformations included in the CARET software [48], first mapped to the human PALS-B12 surface and then to our CIVET template surface. The resulting parcellation is shown in figure 1c, and is also available on request.

Neural activity model and signal expression (a) The neural population model
Our software is flexible concerning the employed model of neural activity. We compute here at each vertex the NPM of Liley et al. [49], as extended by Bojak & Liley [28]. It contains one excitatory (e) and one inhibitory (i) NPM neuron, which represent those populations of real neurons whose coherent activity dominates the macroscopic signal of interest. For example, one can speculate that the EEG is largely due to pyramidal layer V 'output neurons', which form long dendritic bundles [50] that can act as dipolar current sources. The NPM neurons are locally connected to each other (e → i, i → e) and to themselves (e → e, i → i); see figure 3c. NPM self-connections model real neurons of the same type connecting to each other.
The NPM consists of the following ordinary differential equations (ODEs): where l, k = e, i and k lk | e lk =0 ≡ 1 is continuous in e lk . The model parameters could vary from vertex to vertex, but we choose here uniformly those of [51] and set e lk = 0 throughout. 1 Equation (3.1) gives the response of the mean soma membrane potential h k to a sum of post-synaptic potentials (PSPs) I lk . In the absence of input, h k decays exponentially to h r k with characteristic time t k . PSP impact is weighted, with a sign change at the Nernst potentials h lk . The second term, p lk , allows for extracortical input. Here, one could insert thalamocortical loops [52] or specific sensory inputs. We assume simply that extracortical input is noise-like in p ee [28,49]. As a continuous average of neural signals, p ee should be low-pass filtered. We use the spline variant of Catmull & Rom [53] with noise innovations as control points to construct such an input. This allows us to minimize the computational expense for random number generation and spline coefficient computation, which is important since the noise input is computed for every node independently. Catmull-Rom splines are unbiased, interpolating cubic splines with zero local tension and C 1 (first derivative) continuity [54]. The following pseudo-code illustrates how one can interpolate p ee at time steps s = 0, 1, 2, . . ., where time t = s Dt, from the Gaussian white noise rand n sampled every uth time step.

(b) Connectivity and propagation of activity
The third term of equation (3.3), f lk , supports long-range connectivity. We distinguish two kinds; cf. figure 3b. Background connectivity is roughly isotropic, homogeneous and diminishes exponentially with distance [55][56][57], typically leading to waves of cortical activity [58]. This type of connectivity has been commonly used, since it can be approximated by PDEs [38,39] Background synaptic weights w lk,ba sum to one and are multiplied by N a lk,a , the number of synapses formed at vertex a. The firing rates s l,b from vertices b arrive with delays D ba /v lk,ba , where D ba indicates the distance along the cortical surface and v lk,ba the conduction velocity; D ba is estimated as the shortest path through the cortical mesh. We consider here only distances up to a cut-off D c , where exp −L lk,ba D c = 0.1, in order to limit the number of connections to those with significant impact. Long-range connectivity is here considered as exclusively excitatory (N a ik,a ≡ 0), with characteristic decay (1/L ek,ba = 2.5 cm) and cut-off (D c 5.76 cm) distances in the right range for the loss of coherence measured with subdural electrodes [59]. Note that other synaptic footprints can be introduced simply by changing the functional form ofw lk,ba . Specific connectivity is implemented by adding further synaptic weightsŵ ba : where the conduction delays t lk,ba must be given. We will explain our estimation method below. For example, a specific synaptic weightŵ ba = 0.1 means that the added connectivity from vertex b to vertex a has 10% of the strength of the total background connectivity that a is receiving, since b w ba = 1. Specific connectivity can accommodate arbitrary connections between brain regions, and typically will be constructed according to some experimental connectivity matrix. Since the conduction delays D ba /v lk,ba and t lk,ba are not known with great accuracy, we discretize both internally as multiples of the time step Dt = 5 × 10 −5 s. This makes delay bookkeeping much easier. The f lk,a represent averaged and hence continuous signals, for which we again use the spline approach detailed above.
Let us consider f cut = 1.6 kHz as sufficient, then up-sampling u = 5 follows. Thus, we actually need to transmit f lk,a values only every uDt = 0.25 ms. However, the spline needs four control points spaced by three time steps, hence we can only spline delays greater than 3uDt = 0.75 ms. For v ee,ba = 3 m s −1 , the minimum allowed distance is D ba = 2.25 mm. This demonstrates the connection between data transfer and shortest edge length: if we reduce the latter, then we must lower the up-sampling and hence transmit f lk more often.
For parallel computation, we divide up the cortical mesh between compute nodes. To minimize the communication overheads, we compile data into chunks for transmission between nodes. Assume a compute node B has been assigned vertices v b with b ∈ B. A different node A contains vertices v a with a ∈ A. Call B → A ⊂ B the list of vertices for which connection strengths w lk,ba orŵ lk,ba , or both, are non-zero. The compilation of values {s l,b (t)} with b ∈ B → A is what gets transmitted from B to A. To simplify our explanation of how node A distributes the received data, we will assume that data are exchanged at every time step Dt. In reality, only the spline control points are sent. This leads to a significant reduction of traffic, here by a factor u = 5. Every vertex has its own 'delay buffer' for accumulating inputs, with a size set by the maximum (discretized) delay of incoming connections. A pointer indicating 'current time' advances in this buffer and at its end cycles back to the beginning. Distributing {s l,b(t) } now simply reduces to adding N a lk,a bŵ lk,ba s l,b (t) to the position t lk,ba /Dt steps ahead of the 'current time' pointer (modulo the buffer length), and likewise for background connections. When the 'current time' pointer has advanced to a new position in the buffer, the sum of equations (3.6) and (3.7) will have been built up there. The buffer entry is hence added to the system as the current input f lk,a , and then reset to zero, ready to accumulate input again.
The RM does not provide an estimate for fibre length and hence time delays t lk,ba . To calculate an estimate, we assume that vertices are connected by the shortest possible route, which minimizes conduction delays in the brain. Our algorithm first determines all those straight lines between pairs of vertices, which are completely contained within the cortical volume. This yields a network of allowed paths. At this stage any pair of vertices is connected at least via the surface edges, but 'shortcuts' through the volume generally exist. The optimal

(c) EEG and fMRI BOLD signal sources
The details of extracting EEG and fMRI BOLD signals from the NPM simulation have been described in our previous publication; see Bojak et al. [42] for details. Electric dipole strength is here assumed to be roughly proportional to the average excitatory membrane potential h e [61]. A volume conductor model with three compartments was constructed, representing the scalp (conductivity 0.2 S m −1 ), the skull (0.03 S m −1 [62]) and the inside of the skull (0.2 S m −1 ). For piecewise homogeneous volume conductors such as this, the boundary element method can be used to compute the EEG transfer matrix. We assume that fMRI BOLD is driven by glutamate release [63,64]. The neural drive z is hence proportional to the sum f e A ee + f i A ei of excitatory inputs only, cf. equation (3.3). Values f e 0.85 and f i = 1 − f e represent the fraction of excitatory and inhibitory neurons, respectively. We implement a 'balloon-windkessel' model of haemodynamics at every vertex with four ODEs according to Friston et al. [65]. This predicts the local fMRI BOLD contrast y. For the results shown here, a baseline of resting activity is subtracted by hand.

(d) Comparing voxel data with surface simulations
The prediction of voxel values from vertex ones must be informed by anatomy: in the normal direction of the surface we assume here that fMRI BOLD is roughly uniform across the cortical mantle. It is however straightforward to introduce a gradient, e.g. to weigh contributions from the pial surface more strongly. The weighting factors are then determined as follows (cf. figure 5a): if we designate d ij as the vector between a given vertex s i and voxel centre point v j , then n ij is its projection onto the surface normal at s i ; and t ij is a tangential vector uniquely defined by being orthogonal to n ij and coplanar with n ij and d ij . Weights are computed by the Gaussian dispersion 8) where G N,T (x) = exp[−x 2 /(2s 2 N,T )], and we use s N = 2.0 mm and s T = 1.66 mm. The latter equates the Gaussian full width at half-maximum with the average edge length of our cortical surface (3.9 mm). Finally, the projected fMRI BOLD contrast y(v j ) is where n is the number of vertices in the surface mesh. A global normalization restores the original range of fMRI BOLD values. The same approach can be used as well to project fMRI BOLD voxel data onto surface vertices.

Results
The COCOMAC matrix provides basic information about connectivity strength by assigning values 0 (absent), 1 (weak), 2 (moderate) and 3 (strong) based on the reported histochemical staining. Here, we explore the relationship between COCOMAC's anatomical connectivity strength and the effective (causal) connectivity of our model, gauge the strength of specific versus background connectivity, and test the functional visibility of a cortical network consisting of the dorsal (VACd) and ventral (VACv) anterior visual cortex and the frontal eye fields (FEF). As can be seen in figure 1c, these areas form a triangle of connections, all with 'moderate' COCOMAC strength. We also include contralateral connectivity from every vertex to the vertex nearest to its mirror position in the other hemisphere. For the sake of simplicity, we employ only one universal w value for all specific connections. Varying this connection strength parameter will provide a basic test for how significant current experimental and theoretical uncertainties concerning (effective) connectivity are for model predictions of brain dynamics. We have simulated 10 s of brain time forŵ = 0%, 5%, . . . , 100%, i.e. from no specific connectivity to one as strong as the background connectivity. The only input here is the noise shaped as described above, hence all dynamical features at frequencies well below f cut = 75 Hz emerge from the system. The power spectral density of the excitatory mean soma potential h e (simply called 'ECoG' henceforth) was computed for the last 4 s of brain time (2 s Hann window, 50% overlap), avoiding initial transients and averaged over all vertices. The result is shown in figure 6: a strong peak is evident in the delta/theta band, which drifts from 5 to 4 Hz for increasingŵ. There is also a weak alpha peak at 9.5 Hz, which weakens further as it drifts to higher frequencies. Strikingly, upon moving from w = 85% toŵ = 90%, suddenly a strong beta peak appears at around 16 Hz, which forŵ = 100% reaches three-quarters of the maximum power of the slow oscillation. The power of the slow oscillation itself by then has more than doubled. All the described features of the power spectra are due to self-sustaining oscillations, i.e. can be elicited also without the noise drive. We have tested that these results are stable for the simple (forward Euler) integration scheme used here by repeating the calculations with a five times smaller time step.
In figure 7, we compare activity patterns for (a) scalp EEG, (b) ECoG and (c) fMRI BOLD forŵ = 0%, 30%, 60%, 85%, 90%, from top to bottom. In the EEG column, we also see recordings from three 'electrodes' indicated in purple, which show the last second of the simulation. Without specific connectivity, large waves of activity dominate. fMRI BOLD contrast gets sufficiently strong to be detected over experimental noise only between 30 and 60 per cent, and VACd, VACv and FEF also become somewhat visible there. The transition from 85 to 90 per cent appears less dramatic in the activity snapshots than in the power spectrum. However, for 90 per cent and higher we can easily identify the regions of our chosen cortical network, in particular with fMRI BOLD contrast. In summary, our variation ofŵ identified large parameter regions for which brain dynamics would appear qualitatively unchanged in simultaneous EEG/fMRI (30-85% and greater than 90%, respectively), but with the possibility for a transition between these regions induced by a small parameter change (from 85% to 90%).

Conclusions
In the present work and in Bojak et al. [42], we have outlined a complete simulation pipeline for forward predictions of simultaneous EEG and fMRI BOLD. It is based on a well-known NPM formalism [28,49], but uses a cortical tessellation with correct anatomical geometry in a realistic head model. Our software features a method for tracking conduction delays that allows the implementation of arbitrary corticocortical connectivity. We have explained here our pruning procedure for surface meshes, Catmull-Rom splining in noise generation and (multi-parallel) communication, the implementation of specific connectivity and finally the projection of surface fMRI BOLD predictions to voxel data. Furthermore, we have investigated the effects of strengthening specific connectivity. We found that a dynamical regime of dominant slow oscillations that allow the detection of functional networks with fMRI BOLD was qualitatively stable over a large parametric range (specific connectivity strength 30-85% relative to that of the homogeneous and isotropic 'background connectivity'). However, a further small parameter increase (by 5%) would then result in the sudden appearance of additional strong beta oscillations. Consequently, for the time being it appears necessary to adjust connectivity strength carefully with respect to the resulting dynamics, since it is not guaranteed that predictions will be qualitatively similar for two different, but reasonable, parameter choices. We speculate that this difficulty generalizes beyond our current model, and suggest that any dynamical implementation of a connectivity matrix has to be accompanied by a careful exploration of the used connectivity strength.
With respect to frequency content, our simulations may already resemble stage II sleep. However, since our model oscillations are of purely cortical origin, lacking both realistic thalamic input and sleep cycle modulation, the current setup provides no direct insight into the mechanisms of sleep. Importantly, however, the NPM and connectivity can be easily modified or exchanged. Hence the sleep NPMs mentioned in §1 (cf. also Robinson et al. [66] in this Theme Issue) could be integrated with our software to allow realistic predictions of the activity of the sleeping brain as observed with (simultaneous) EEG and fMRI BOLD.