simulator.tex

The simulator in TVB resembles popular neural network simulators in 
many fundamental ways, both mathematically and in terms of informatics 
structures, however we have found it necessary to introduce auxiliary
concepts particularly useful in the modeling of large scale brain 
networks. In the following, we will highlight some of the interesting
principles and capabilities of TVB's simulator and give rough characterization
of the execution time and memory required in typical simulations.

\subsection{Node dynamics}

	In TVB, nodes are not considered to be abstract neurons nor necessarily
	small groups thereof, but rather large populations of neurons. Concretely,
	the main assumption of the neural mass modeling approach in TVB is that
	large pools of neurons on the millimeter scale are strongly approximated
	by population level equations describing the major statistical modes of
	neural dynamics \citep{Freeman_1975book}. Often, averaging techniques are
	employed, though techniques retaining several modes have been developed
	\citep{Stefanescu_2008, Stefanescu_2011}. Such an approach is certainly not
	new; one of the early examples of this approach consist of the well known
	Wilson-Cowan equations \citep{Wilson_1973}. Nevertheless, there are
	important differences in the the assumptions and goals from modeling of
	individual neurons, where the goal may be to reproduce correct spike
	timing or predict the effect of  a specific neurotransmitter. A second
	difference lies in coupling: chemical coupling is often assumed to be
	pulsatile, or discrete, between neurons, whereas for mesoscopic models it is considered
	continuous. Typically the goal of neural mass modeling is to study the
	dynamics that emerge from the interaction of two or more neural masses and
	the network conditions required for stability of a particular
	spatiotemporal pattern. In the following, we shall  briefly discuss some
	of the models available in TVB.

	As we have noted, many neural mass models have been developed. One of
	the more prominent examples in the systems neuroscience literature is 
	the Jansen-Rit model of rhythms and evoked responses arising from
	coupled cortical column \citep{Zetterberg_1978, Jansen_1995,David_2004, Spiegler_2010}. 
	Advantages of the Jansen-Rit model stem from the connection made
	between empirical studies of neural tissue and the model's parameters, 
	making it easier in certain cases to make concrete predictions about
	the relation between a dynamical regime and its neurobiological 
	mechanism. However, because the form of the model used often employs
	at least six dimensions, it is not always clear how it should be analyzed or
	visualized. Lastly, the model requires frequent computation of exponentials,
	requiring considerable computational time. 

	For these reasons, it is often desirable to have a simpler mathematical 
	model, which can reproduce the same qualitative phenomena as other 
	models, implemented with fewer and simpler equations. Such is the motivation
	for the generic two-dimensional oscillator model 
    provided by TVB 
    \citep{strogatz2001nonlinear, guckenheimer1983nonlinear}.
	This model can produce oscillations which are damped, spike-like or 
	sinusoidal in nature. While these alone are not interesting, they 
	permit the study of network phenomena, such as synchronization of rhythms
	or propagation of evoked potentials, while requiring less time to simulate.

	However, the modeler's goals may not lead to either the Jansen-Rit
	 or generic 2D oscillator, so several other mass models are
	provided by TVB: the previously mentioned Wilson-Cowan description of
	functional dynamics of neural tissue \citep{Wilson_1972}, the Kuramoto
	model describing synchronization \citep{Kuramoto_1975, Cabral_2011}, two
	and three dimensional mode-level models describing populations with
	excitability distributions \citep{Stefanescu_2011, Stefanescu_2008}, a
	reduction of \citeauthor{Wong_2006}'s (\citeyear{Wong_2006}) model as presented 
    by
	\cite{Deco_2013} and a lumped version of 
    Liley's model \citep{Liley_1999, Steyn-Ross_1999} 
    are among the available models in TVB.

	Again, should any of these be insufficient, a new model can be implemented
	with minimal effort by subclassing a base \texttt{Model} class and
	providing a  \texttt{dfun} method to compute the right hand sides of the
	differential  equations. Please refer to the \texttt{tvb.simulator.models}
    module of the main scientific library or the \texttt{contrib} folder
    for
	examples. In the \texttt{contrib} folder, models from the work of 
	\cite{Larter_1999, Breakspear_2003, Morris_1981, Hindmarsh_1984, Brunel_2001}
	have been implemented.

\subsection{Network structure}

	The network of neural masses in TVB simulations directly follows from  a
	pair of geometrical constraints on cortical dynamics. The first is the
	large-scale white matter fibers that form a non-local and heterogeneous
	(translation variant) connectivity, either measured by anatomical tracing
	(CoCoMac, \cite{Koetter_2004}) or diffusion-weighted imaging 
    \citep{Hagmann_2008, Honey_2009, Bastiani_2012}.
    The second is that of horizontal projections
	along the surface, which are often modeled with a translation invariant
    connectivity kernel, approximating a neural field; though as with other
    parameters in the simulator, spatial inhomogeneity is supported as well.

	\subsubsection{Large-scale connectivity}

	The large-scale region level connectivity at the scale of centimeters,
	resembles more a traditional neural network than a neural field, in that,
	neural space is discrete,  each node corresponding to a neuroanatomical
	region of interest, such as V1, etc. It is at this level that inter-regional 
	time delays play a large role, whereas the time delays due to 
	lateral, local projections are subsumed under the dynamics of the node.

	It is often seen in the literature that the inter-node coupling functions
	\textit{are} part of the node model itself. In TVB, we have instead 
	chosen to factor such models into the intrinsic neural mass dynamics, where each 
	neural mass's equations specify how connectivity contributes to the
	node dynamics, and the coupling function, which specifies how the activity
	from each region is mapped through the connectivity matrix. Common coupling 
	functions are provided such as the linear, difference and periodic functions
	often used in the literature.

	\subsubsection{Local connectivity}

	The local connectivity of the cortex at the scale of millimeters provides
	a continuous 2D surface along horizontal projections connect 
	cortical columns. Such a structure has previously been modeled by
	neural fields \citep{Amari_1977, Jirsa_1997, Liley_1999}. In TVB, a cortical mesh, 
	as obtained from structural MRI data and simplified, provides a spatial 
	discretization on which neural masses are placed and connected with a
	local connectivity kernel, itself only a function of the geodesic distance
	between the two masses. This is considered to provide a reasonable
	approximation of a neural field, the appropriateness of which depends on the properties
	of the mesh and the imaging modalities that sample the activity simulated
	on the mesh \citep{Spiegler_2013}. In fact,
	the implementation of the local connectivity kernel is such
	that is can be re-purposed as a discrete Laplace-Beltrami operator,
	allowing for the implementation of true neural field models that 
	use a second-order spatial derivative as their explicit spatial term.

	TVB currently provides several connectivity kernels, of which a Gaussian
	is one commonly used. Once a cortical surface mesh 
	and connectivity kernel and its parameters are chosen, the geodesic
	distance (i.e. the distance along the cortical surface) is evaluated
	between all neural masses \citep{Mitchell1987}, and a cutoff is chosen
	past which the kernel falls to 0. This results in a sparse matrix that 
	is used during integration to implement the approximate neural field. 

\subsection{Integration of stochastic delay differential equations}

	In order to obtain numerical approximations of the network model 
	described above, TVB provides both deterministic and stochastic
	Euler and Heun integrators,
	following recent literature on numerical solutions to stochastic
	differential equations \citep{Kloeden_1995,Mannella_2002,Mannella_1989}.

	While the literature on numerical treatment of delayed or 
	stochastic systems exists, it is less well known how to treat 
	the presence of both. For the moment, the methods implemented by TVB
	treat stochastic integration separately from delays. 
	This separation coincides with a modeling assumption that in
	TVB the dynamical phenomena to be studied are largely determined
	by the interaction of the network structure and neural mass dynamics, 
	and that stochastic fluctuations do not fundamentally reorganize the
	solutions of the system \citep{Ghosh_2008,Deco_2009,Deco_2011,Deco_Senden_2012}.

	Due to such a separation, the implementation of delays in the
	regional coupling is performed outside the integration step,
	by indexing a circular buffer containing the recent simulation 
	history, and providing a matrix of delayed state data to the 
	network of neural masses. While the number of pairwise
	connections rises with $n_{region}^2$, where $n_{region}$ is
	the number of regions in the large-scale connectivity, 
	a single buffer is used, with a shape
	$(horizon, n_{cvar}, n_{region})$ where $horizon = max(delay) + 1$,
	and
	$n_{cvar}$ is the number of coupling variables. Such a scheme helps 
	lower the memory requirements of integrating the delay equations.

\subsection{Forward solutions}

	A primary goal of TVB is not only to model neural activity itself
	but just as importantly the imaging modalities common in human 
	neurosciences, using so-called forward solutions, which allow for
	the projection of neural activity into sensor space. To account
	parsimoniously for other ways in which simulated data might be saved, 
	such as simple temporal averaging, we refer to each of these simply as 
	\textit{Monitors}, which take as input neural activity and 
	output a particular projection thereof. In most cases, this 
	takes the discrete-time form of

	\[ \hat{y}[j, t] = \sum_{i=1, \tau=1}^{N_W, N_k} W[j, i] K[\tau] y[i, t-\tau] \]

	\noindent where $y[i, t]$ is the amplitude of the $i^{th}$ neural mass at time
	$t$, $K[\tau]$ is a temporal kernel, and $W[j, i]$ is a spatial kernel,
	usually projecting the state variable of interest of the $i^{th}$ 
	neural mass to the $j^{th}$ sensor. 

	Where necessary for computational reasons, monitors employ more than 
	one internal buffer. The fMRI monitor is one 
	example: given a typical sampling frequency of simulation may be upward of 
	64 kHz, and the haemodynamic response function may last several seconds, 
	using only a single buffer could require many gigabytes of memory for the 
	fMRI monitor alone. Given that 
	the time-scale of simulation and fMRI differ by several orders of magnitude, 
	the subsequent averaging and downsampling is justified. 

	In the cases of the EEG and MEG monitors, $K$ implements a simple
	temporal average, and $W$ consists of a so-called lead-field matrix as typically
	derived from a combination of structural imaging data of the patient, 
	which provides the locations and orientations of the neural sources, and the locations
	and orientations of the EEG electrodes and MEG gradiometers and magnetometers. 
	As the development and implementation of such lead-fields is well developed
	elsewhere \citep{Sarvas_1987,Hamalainen_1989,Jirsa_2002,Nolte2003,Gramfort_2010}, TVB provides access
	to the well-known OpenMEEG package, however, the user is free to provide 
	their own.

\subsection{Performance}

	A primary goal of the simulator is to be available as a pure Python package,
	and secondarily, to be fast enough. We have not found it useful to 
	develop theoretical estimates of the time and space complexity of the 
	algorithms, given that much of the heavy lifting is already done in native
	code by NumPy and other standard libraries. Instead, 
	in the following, we profile a set of eight characteristic simulations
	on both memory use, specifically the heap size as measured by Valgrind's 
	\texttt{massif} tool \citep{nethercote2007valgrind}, 
    and function timing as measured by the 
	\texttt{cProfile} module of the standard library. 
	
	Measurements were
	performed on an HP Z420 workstation, with a single Xeon E5-1650
	six-core CPU running at 3.20 GHz, L1-3 cache sizes 384 KB, 1536 KB
	and 12 MB respectively, with main memory 4 x 4 GB DDR3 at 1600 Mhz,
	running Debian 7.0, with Linux kernel version 3.2.0-4-amd64. 
	The 64-bit Anaconda Python distribution was used with additional Accelerate
	package which provides acceleration of common routines based on the 
	Intel Math Kernel Library. A Git checkout of the trunk branch of TVB 
	was used with SHA \texttt{6c644ab3b5}.

	Eight different simulations were performed corresponding to the combinations of
	either the generic 2D oscillator or Jansen-Rit model, region-only
	or use of cortical surface, and two conduction speeds, $v_c = 2.0$ and
	$v_c = 20.0$ (m/s). In each case, a temporal average monitors at 512 Hz
	is used, and the results are discarded. The region-only simulation was
	run for one second while the surface simulation was run for 100 ms. 

	\note[mw]{Table of profiling results to go here. Profiling has been done, 
	table will be added soon. Nothing surprising here.}

%\subsection{Acceleration with C \& CUDA}

%	\input{c-cuda.tex}