  This paper introduces methods and a novel toolbox that efficiently integrates
any high-dimensional Neural Mass Models (NMMs) specified by two essential
components. The first is the set of nonlinear Random Differential Equations of
the dynamics of each neural mass. The second is the highly sparse
three-dimensional Connectome Tensor (CT) that encodes the strength of the
connections and the delays of information transfer along the axons of each
connection. Semi-analytical integration of the RDE is done with the Local
Linearization scheme for each neural mass model, which is the only scheme
guaranteeing dynamical fidelity to the original continuous-time nonlinear
dynamic. It also seamlessly allows modeling distributed delays CT with any
level of complexity or realism, as shown by the Moore-Penrose diagram of the
algorithm. This ease of implementation includes models with distributed-delay
CTs. We achieve high computational efficiency by using a tensor representation
of the model that leverages semi-analytic expressions to integrate the Random
Differential Equations (RDEs) underlying the NMM. We discretized the state
equation with Local Linearization via an algebraic formulation. This approach
increases numerical integration speed and efficiency, a crucial aspect of
large-scale NMM simulations. To illustrate the usefulness of the toolbox, we
simulate both a single Zetterberg-Jansen-Rit (ZJR) cortical column and an
interconnected population of such columns. These examples illustrate the
consequence of modifying the CT in these models, especially by introducing
distributed delays. We provide an open-source Matlab live script for the
toolbox.
