Computation of the MEG response from the model
Dynamical and topographical parameters
To this end, we have focused on the description of the system dynamics of the AC from which the decay constants, angular frequencies, and damping frequencies emerge. However, the MEG response actually originates from a mixture of a dynamical and a topographic part. The dynamical element is captured by a set of system parameters \(\bf{P}\) which are essentially determined by the various weight matrices, see Equations (\ref{eq:udsol})–(\ref{eq:freq}). The topographical element reflects a non-dynamical modulation of the signal which accounts for the topography of the primary currents and is condensed into a matrix \(\bf{K}\). It includes the effects due to the type of connection represented by each synapse, i.e. whether it is a feedforward or feedback connection, the distance between the sensor and the cortical field in which the respective connection is located, as well as the subject-specificity of the topography of the cortical surface.
The \(\bf{K}\) matrix displayed in Figure \ref{fig:kmat} describes the sign and magnitude of the contributions of each synaptic connection to the MEG signal. The underlying assumption is that area-to-area feedforward signals arriving in the middle layers (predominantly in layer IV) target proximal locations of the pyramidal dendrites and result in a current flow which points upward; this is represented by light green elements with negative values in Figure \ref{fig:kmat}. Area-to-area feedback signals target the upper layers (predominantly layers I and II) and therefore arrive at apical dendritic locations resulting in a current flow downward, indicated by purple elements with positive values in Figure \ref{fig:kmat}. Thus, feedforward and feedback signals lead to primary currents pointing in opposite directions, towards and away from the cortical surface, respectively \citep{Ahlfors2015}. This, in turn, means that the topography of the cortical surface affects these signals through the distance between primary currents and the sensor, and the direction of the primary current with respect to the surface of the head \citep{Hamalainen1993}. We point out that arriving at a “true” \(\bf{K}\) matrix would involve a complete description of the cortical surface of the individual subject, a complete anatomical layout of the cortical fields, as well as a complete description of the primary currents at each cortical location. This is practically equivalent to having access to the inverse solution via a complete anatomical picture of the subject-specific AC. This is clearly an impractical and unrealistic requirement. Thus, we lowered our expectations and aim for insights about the dynamics of the system without having to solve the inverse problem.
Simulations of MEG responses
The strength of the excitatory connection between two columns can either enhance (\(W_{\tiny\mbox{ee}}\)) or weaken (\(W_{\tiny\mbox{ie}}\)) the MEG response; these two effects are intertwined in the matrix \(W_{\tiny\mbox{ee,lat}}\). The computation of the MEG response is based on the input currents which can be treated as the weighted sum of spiking rates which, in the analytical solution, equals the weighted sum of the excitatory state variables: \(\sum_{ij}W_{{\mbox{\tiny{ee,lat}}}_{ij}}u_{j}(t)\). Figure \ref{fig:MEG_wave}A shows the respective waveform for the example of underdamped normal modes and coupled state variables presented in Figure \ref{fig:normod-waves}. Taking into account the topographical information embedded in the \(\bf{K}\) matrix, the MEG response \(R(t)\) is computed as the sum of the weighted excitatory (postsynaptic) input currents multiplied with the \(\bf{K}\) matrix \citep{inproceedingsMay2002}:
\begin{equation}
\label{eq:MEG}
\label{eq:MEG}R(t)=\sum_{ij}K_{ij}W_{{\mbox{\tiny{ee,lat}}}_{ij}}u_{j}(t).\par
\\
\end{equation}
The thick blue line in Figure \ref{fig:MEG_wave}B shows the corresponding MEG waveform computed with Equation (\ref{eq:MEG}). Compared to the waveform of the sum of input currents, the overall morphology of the waveform is largely preserved. The MEG response, however, includes now the P50, a small negative deflection prior to and clearly discernible from the M100. Further, the M100-peak latency is shifted to larger values.
We also performed simulations to verify the effect of different input frequencies (associated with different input columns of the three core fields) on the MEG response. The corresponding results are also shown in Figure \ref{fig:MEG_wave}B. Here, in the first simulation, the input signal from thalamus was simultaneously targeting the first column of the three core fields, the second column in the second simulation and so on, until all 16 columns have been separately targeted by the input. We observed a remarkable stability of the morphology among the resulting 16 waveforms, with a distinct dependence of the M100-peak magnitude, and less so of the peak latency, on the input column (i.e. input frequency). This dependence of the M100-peak latency on the input frequency is displayed in the inset of Figure \ref{fig:MEG_wave}B. Albeit minimum and maximum peak latencies differ only by a few milliseconds, a U-shaped behaviour can be observed which is in qualitative agreement with experimental observations reported, for example, by \citet{Roberts1996} and \citet{Stufflebeam1998}.
The P50 deflection emerged in simulations of the MEG response only when the topographical information comprised in the \(\bf{K}\) matrix was taken into account. In fact, the \(\bf{K}\) matrix allows investigating generic effects of feedforward and feedback connections on the overall MEG waveform, in particular on the various positive and negative deflections. In Figure \ref{fig:fb_ff}A, the strength of feedforward connections is stepwise altered from \(-5\) (thick blue line) to \(-50\) (red line) while keeping the weights of feedback and intra-field connections constant.
We found systematic changes of the peak magnitude and latency of both the P50 and the M100, albeit with a much more pronounced amplitude effect on the P50. By contrast, varying the strength of the feedback activity in the same absolute range, from 5 to 50, while keeping the weights of the feedforward and intra-field connections constant
results in a monotonous increase of the M100-peak magnitude with only minor changes of the P50, see Figure \ref{fig:fb_ff}B.
Information on the loci of MEG response generation can be derived from simulated MEG waveforms by splitting them up into contributions from the underlying AC fields and areas. This is shown in Figure \ref{fig:MEG_field_area} where the MEG waveform (thick blue line) displayed in Figure \ref{fig:MEG_wave}B is broken up into the individual waveforms of the 13 fields (Figure \ref{fig:MEG_field_area}A) and into the summed responses for the three AC areas, i.e. for core, belt and parabelt (Figure \ref{fig:MEG_field_area}B). We found distinct area-specific field contributions to the overall MEG response, with decreasing M100 peak magnitudes and increasing peak latencies from core to belt to parabelt. We also note that the generation of event-related responses is accounted for in terms of the internal structure of auditory cortex and its dynamics. In this account, each peak and trough of the event-related response (e.g., the M50, M100, M200) is not due to a dedicated response generator but, rather, arises out of the network properties of the entire AC.
This notion is supported by further simulations where we computed response waveforms of the individual AC areas as well as of the overall MEG response using different connection patterns and weight values for the two parabelt fields while keeping all other weights of \(\widetilde{W}_{\tiny\mbox{ee,lat}}\) constant. On the supposition that the MEG waveform reflects network properties of the entire AC rather than the superposition of distinct focal response generators, we anticipated that even changes of parabelt connections only, i.e. alterations of weight values at the fringe of AC, should affect the MEG response. Figure \ref{fig:pb_change_MEG} is a compelling proof-of-concept for this idea. It is an example of the impact of different inter- and intra-field connection patterns of the parabelt fields on the response waveforms of all areas and the overall MEG. Figure \ref{fig:pb_change_MEG}A represents the connections of the parabelt fields originally introduced in Figure \ref{fig:AC_w_ee}B; Figure \ref{fig:pb_change_MEG}B is a modification of it. Note that the general connection pattern is preserved by and large, and only the parabelt weights have been altered; all other weights of \(\widetilde{W}_{\tiny\mbox{ee,lat}}\) were unchanged. Figures \ref{fig:pb_change_MEG}C, D, and E show the resulting areal waveforms for the original (solid lines) and the modified (dashed lines) case for the core, belt, and parabelt, respectively. In all three areas, the mere modification of the parabelt connections entails significant changes in the respective waveforms, resulting, for example, even in a larger M100-peak magnitude of the belt waveform compared to the one from the core, or an opposite behavior of the two waveforms in the parabelt. The summed MEG response of the three areas is shown in Figures \ref{fig:pb_change_MEG}F. In the example considered here, the modification of the parabelt connections strongly affects the MEG response and leads to a much broader waveform with a larger M100-peak amplitude and latency and a vanishing of the P200 peak, i.e. a flattening of the waveform beyond the M100 deflection.
Separating dynamics from topography
In further simulations we tested whether the dynamical and topographical parameters differently affect the MEG response, i.e. whether the MEG signal can be decomposed in such a manner that the dynamical system can be identified and separated from the effect of K. The overarching question here is whether there are aspects of the MEG signal which depend only on the dynamics of the system but not on the K modulation. Specifically, we examined the MEG response under two opposing conditions: 1) The dynamical parameters \(\bf P\) were kept constant and the values of the feedforward, feedback and intra-field weights included in the \(\bf K\) matrix were altered (Figs \ref{fig:Rand_P_K}A-D), and 2) the elements of the \(\bf K\) matrix were kept constant, but the weight matrices forming \(\bf P\) were modified (Figs \ref{fig:Rand_P_K}E-H). Figure \ref{fig:Rand_P_K} clearly shows that the randomization of the \(\bf P\) parameters leads to a much larger variation of the MEG response compared to the randomization of \(\bf K\). This does not only hold for the morphology of the MEG signal (compare Figs \ref{fig:Rand_P_K}A, E), but is also visible in the respective spectral analysis performed by means of FFT (Figs \ref{fig:Rand_P_K}B, F), and in the attenuation constant derived from the waveforms (Figs \ref{fig:Rand_P_K}C, G and Figs \ref{fig:Rand_P_K}D, H).
Note that the separability of the effects of \(\bf P\) and \(\bf K\) is a simplification stemming from the initial approximation of describing each column with only two state variables (\(u\) and \(v\)). It assumes that each column has a 2-D state which evolves independently of the locus of the input, and that this locus only affects the generation of the MEG signal. In reality, each column is a complex dynamical system of its own with multiple state variables, layers, and input loci, and where the input location probably affects the evolution of the signal (i.e., the impulse response function depends on impulse location). How seriously our approximation deviates from a biophysically realistic description of the column would involve extensive numerical studies beyond the scope of the current efforts. Notwithstanding that our formulation of \(\bf K\) is (necessarily) coarse, our simulations demonstrate that we have access to a dynamical-systems description of AC which bypasses \(\bf K\). This means the MEG signal can be decomposed in such a manner that the dynamical system can be identified and separated from the effect of \(\bf K\). This is of particular importance for our future approaches of fitting the analytical solutions to experimental MEG data to derive subject-specific model parameters. These findings reinforce our assumption that this fitting procedure captures the dynamical descriptors (normal modes) even when the true \(\bf K\) matrix remains unknown.