Accurate and efficient Simulation of very high-dimensional Neural Mass Models with distributed-delay Connectome Tensors

: 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. To date, simplistic assumptions prevail about delays in the CT, often assumed to be Dirac-delta functions. In reality, delays are distributed due to heterogeneous conduction velocities of the axons connecting neural masses. These distributed-delay CTs are challenging to model. Our approach implements these models by leveraging several innovations. 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. This semi-analytic LL integration is highly computationally efficient. In addition to this, a tensor representation of the CT facilitates parallel computation. 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.


Introduction
Neural Mass Models (NMM) are essential for exploring mesoscopic nonlinear brain dynamics.Obtained as meanfield approximations of ensembles of neurons, they retain sufficient simplicity to allow reasonably efficient implementation and analysis.Close cousins of these models are Neural Field models, which are not dealt with in this work but can be approximated by a sufficiently high dimensional NMM.For a review of both types of models, see (1).
NMMs serve a different purpose than models targeting the microscopic level.Programs such as GENESIS (2,3) and NEURON (4,5) allow in-silico modeling realistic neuronal aggregates.An outstanding example of this line of work is the Blue Brain (BB) (6)(7)(8) project, which aims to construct a biophysically detailed simulation of a cortical column with the most exquisite anatomical definition.This microscopic level of detail needs to be circumvented in particular circumstances, such as comparing brain models with mesoscopic neuroimaging data.The reasoning here is that a coarser level of modeling retains the essential emergent properties of neural aggregates but allows more straightforward calculations, especially for solving inverse problems.There are excellent implementations of NMMs, examples being Dynamic Causal Modelling (DCM) (9,10), The Brain Dynamics Toolbox (11), The Virtual Brain (TVB) (12)(13)(14)(15), among others.This paper introduces methods and a toolbox for efficiently integrating any high-dimensional Neural Mass Model (NMM).To understand our paper's motivation, consider specific modeling and algorithmic issues that currently limit NMM capabilities.Two essential components specify these NMM.
• The dynamics of each Neural Mass: The first component is the set of nonlinear Deterministic, Stochastic, or Random Differential Equations (RDE) that specify the dynamics of each neural mass (or node).• The second component is the highly sparse three-dimensional Connectome Tensor (CT).This tensor encodes the strength and the delay of information transfer for each connection between nodes.The three modes of the CT are emitter nodes, receiver nodes, and the distribution of delays between nodes.
Figure 1 shows the general structure of the Connectome Tensor 1 .Most work on functional connectomics collapses the CT across the temporal domain, reducing it to a connectivity matrix.Such a matrix formulation is adequate for anatomical connectivity but ignores temporal delays of signal transmission along connection pathways.The low temporal resolution of fMRI has perpetrated such an approach but is undoubtedly inadequate to disentangle causal relationships among brain regions (16).
The complexity of both node dynamics and the CT have generated conceptual and algorithmic.Here we address specific algorithmic issues that have been partially solved in successive stages.We now summarize these and explain the contributions of this paper: whose solution might enhance current NMM software capabilities.These concerns and the solutions explored in this paper are: Lack of Fidelity to the original dynamical system: NMM integrates nonlinear deterministic (DE), stochastic, or random differential equations (SDE or RDE, respectively).Some current NMM implementations use standard solutions of these equations, sometimes blithely using deterministic techniques, which gives erroneous results.In the best cases, methods for SDE are employed, such as Euler-Maruyama, Heun, and 4 th order Runge-Kutta methods (17,18).These methods do not entirely preserve the original dynamics of the continuous system.Also, standard solutions become calculations frequently explosive (19)(20)(21) with probability 1.Also, for stochastic systems and small values of the integration step, the resulting Markov chain is non-stationary for nonlinear functions (22,23).A solution to these problems was introduced in (24), employing the Local Linearization Method (LLM) for the solution of SDE or RDE  (20,25,26).This method guarantees both fidelity and stability to obtain a stationary Markov chain by interaction, which guarantees stationary Markov chains and non-explosive discrete-time dynamical systems, thus providing discrete-time mimicking of the original continuous-time dynamic system (27).LLM was applied for the first time to integrate NMM and parameter estimation of a neural state-space model in (28).Since then, it has been adopted in DCM (29).

Excessive computational complexity:
The LL integration technique has a very high computational cost.Each integration step requires the calculation of matrix exponentials whose dimensionality is of the order of the number of neural masses.Despite the emergence of optimized methods for this calculation based on Krylov space projection (30), the numerical burden was still excessive.It limited the use of LL procedures to small-sized neural networks (20).In the paper (31), computational complexity was tackled by taking advantage of the transmission delays between neural masses.It is therefore essential to decrease the integration step of the NMMs.This increased sampling rate also increases computational cost that grows as higher-dimensional models are explored but compounded when solving inverse neuroimaging problems.Inverse problems based on NNM dynamics require iterative recalculation of state-space trajectories.
Fortunately, with a sufficiently small integration step, the input to each neural mass is already known.Thus, the solution of the RDE can be carried out for each mass in isolation.In this case, the integration for each neural mass can be solved symbolically and expressed explicitly.This tactic allows orders of magnitude acceleration of computation.The input to each neural mass can be expressed as linear combinations of time-shifted versions of the output of the other neural masses, weighted by the appropriate connectivities.This approach was used to simulate realistic, largescale EEG and fMRI realizations in real-time (32).Here the inter-column connectivities of each Zetterberg-Jansen-Rit (ZJR) cortical column were fixed, and the long-range inter-areal connectivities were obtained from Diffusion-Weighted Imaging (DWI) (32).Despite its computational advantages, this software has not been frequently used due to its complexity and lack of detailed algorithm descriptions.It was also based on using a single conduction velocity of axonal communication between nodes.

Difficulty in dealing with realistic profiles of axonal time delays:
When considering the delays in transmission information from one node to another, it is well known that the extreme diversity of myelination (or lack of) and axonal fiber diameters lead to a distribution of conduction velocities (and therefore delays) (33)(34)(35)(36)(37).The complexity of solving delay differential equations has led many teams to ignore them, assuming instantaneous axonal propagation, thus leading to systems of simultaneous differential equations.DCM finesses the issue of delays by tampering with the length of the modeled postsynaptic potential (38)(39)(40).Other groups, including TVB, pose axonal delays as Dirac-delta to facilitate solutions (41,42).This Dirac-delta delay was also used in (32), where a single conduction velocity is associated with each connection.
What we overlooked in our earlier formulation was that incorporating distributed delays could be easily achieved within the RDE-algebraic approach described in (32).We intend to remedy this omission in this paper.As mentioned previously, with a sufficiently small integration step, the inputs to each neural mass may be treated as already available, paving the way for totally general distributed transmission delays.Thus, instead of specifying connectivities between nodes as having a fixed delay, we reformulate the discrete integration scheme as a tensor operation with the three-dimensional sparse CT specifying the emitter neural masses projecting to receiver masses with unrestricted time lag distribution.

Numerical efficiency:
As we have already stated, it is essential to decrease the integration step for the neural mass models to use the semi-analytic integration.Thus, some gains in efficiency are now offset by an increased computational cost due to more integrations.This burden grows even more as higher-dimensional models are explored.The problem is compounded when solving neuroimaging inverse problems based on NNM dynamics requiring iterative recalculation of state-space trajectories.
In this paper, we present theory and software to extend the random differential-algebraic formulation of NMM (already formulated (32)).Here, we perform in tandem, with considerable computational advantage, i) Symbolical integration of the LLM integration step for each separate neural mass and ii) Calculation of inputs to each neural mass with the highly parallelizable CT operations.We make the Matlab Live Script for the algorithm available.Additionally, we illustrate the use of the method in examples that, without pretension to neurobiological fidelity, highlight the use of the algorithm.
We show that by decreasing the integration step, inputs to each neural mass may be treated as already available, thus paving the way for totally general distributed transmission delays.We can now reformulate the discrete integration scheme as a tensor operation with the three-dimensional sparse connectivity of emitter neural masses projecting to receivers with any given lag distribution.

Model Description
Since our formulation depends on tensor operations, we briefly describe our notation summarized in Tables 1 and  2.

Product Description
Hadamard product of Tensor and matrix   Components of the state vector ( ) The nonlinear sigmoid function that converts postsynaptic potentials to spiking rate The maximum amplitude of the PSP Observation equation: ( ) ( ) The state vector of all the neural masses is ( ) which determines the observations ( ) x at discrete time instants with t ξ being the measurement noise.The parameters  determine the properties of the neural masses and the dynamical noise driving the system ( ) t  .Some definitions: • We partition the state vector variables for each neural mass as the set of receiver nodes.
• The set of all neural masses is We can now define the connectivity matrix . It encodes the synaptic strength of the connections between the emitter and transmitter node.This matrix is the usual connectivity matrix of many models (17,31,(44)(45)(46)(47)(48)(49) that do not explicitly model neural transmission delays.
The primary purpose of this paper is to incorporate transmission delays into our models.Towards this end, we define the interval The D tensors used in this paper are:  .In this case, e r d is the length of the fiber system connecting nodes e x and r x a single axon propagation speed.Examples of this approach can be found in (31,50).This approach ignores the diversity of axon diameter and degree of myelination in reality but is simple to estimate.
• A more realistic , where  is an integration constant to ensure that ( )

3
; er p u d is a density, , and e r d is defined as above.This distribution function was proposed by Nunez (51) based on an experimentally determined distribution of corticocortical axon diameters (51).
Other delay models are possible and easily incorporated into our framework.
We are now ready to state our Neural Mass Model with distributed-delay Connectome Tensor expressed as where − x is the delayed system vector and K is the CT defined above.
We can now explain the most significant simplification underlying our proposed algorithms, and it is that we can perform many computations for each neural mass separately.To proceed, we need the following definitions: the M T  matrix of delayed state vectors for each neural mass where x is the input to a neural masses j of the past activity of all other neural masses.
With the definitions in place, we can define our simplified neural mass model as This formulation substitutes the set of random differential delay equations (4) with a set of random differentialalgebraic equations.This approach generalizes the differential-algebraic interpretation from (32), introducing generalized distributed delays and making calculations more accessible and less time-consuming.

Numerical Integration
To deal with numerical computations, we must discretize the system of equations (5).We discretize the connectivity tensor by defining an integration step t  .T is now sampled at multiples of t  .This discretization allows us to redefine  ( )  is the integration interval.
To integrate equation ( 5)b) we apply the Local Linearization Method (LLM) (24,32,52,53), which overcomes the limitations of the classical methods mentioned in the introduction.This method is applicable regardless of the complexity of the neural mass in question.Our approach constructs an equivalent locally linear system for the interval ( )  as described in (26).
The resulting discretization ( ) ( ) of the model is given by The integration starts with the initial value ( ) We obtain the integration function  is theoretically based on the use of the matrix exponential (54,55).However, we efficiently carry out this computation through a specific block triangular matrix combining various specifically designed submatrices.See (56) for the method description.Consequently, the discrete state-space model is simply where 0 ,1 Henceforth, we specialize the general model in equations (5) to the case where each neural mass follows the dynamics stated in the Zetterberg (45)-Jansen-Rit (57) (ZJR) NMM family of models.However, we note that our proposals may be used with other formulations.We note that there are two standard formulations for ZJR NMM.One uses kernels as done in (58)(59)(60)(61)(62)(63).The other form uses the model's state-space form (32), which we follow in this work.

Model for a single neural mass and algebraic integration
For a single ZJR neural mass model formulation, equation ( 5)b) can be rewritten as a first-order system of two differential equations with the vector input ( )  (64).Therefore, following the statement of equation ( 5)b,c), the model for a single neural mass is formulated as where ( ) j y t is the primary variable, which represents the PSP output of a single neural mass and ( ) t is its derivative.
The system parameter vector for a single neural mass is ( ) , , , The values of these parameters determine the system response of the neural mass (See Table 4).The system has two inputs: the noise process ( ) and an external input process ( ) j y t (from other neural masses).
We show the block diagram of a single neural mass in Figure 1.
where 0 e is the maximum average firing rate of the neural populations, 0  is the PSP corresponding to 0 e , and a is the slop of the sigmoid.The values of these parameters depend on whether the PSP is excitatory (EPSP) or inhibitory (IPSP) and can be found in Table 3 (66-68).
Since the system is small enough to be solved symbolically, we carried out the algebraic derivation of the LL integration step explicitly with a Live Script using the MATLAB Symbolic Toolbox.The output of these operations is shown in the expressions (12)-Error!Reference source not found.,where we group the random system input into the variable ( ) Furthermore, the full Jacobian of the system is determined substituting the expressions (5)b), (12) and Error!Reference source not found.obtained by algebraic derivation as follows In the next step of the LLM, after applying the exponential matrix, is a regrouping of terms was made, which simplifies and facilitates the compression and manipulation of the following numerical expression ( ) ( ) ( ) ( ) where j A and j B are constant matrices with the PSP factors for an arbitrary neural mass j M  , I is the identity matrix, and ( ) j t e contains the random inputs to the system and the sigmoidal function (32).
Note that the single neural mass can be considered from now on as a "black box", a node of any ZJR-based neural mass model which can ignore the internal details (such as the mechanics of the numerical integration).It remains now to assemble different models with these nodes.

Single Zetteberg-Jansen-Rit (ZJR) Cortical Column Level
Our model framework now can describe the original Zetterberg -Jansen-Rit (ZJR) model (46,47) for a single cortical column.The ZJR column is a system composed of three neural masses (Pyramidal: ).Each of these is a first-order system of two differential equations ( 9), thus resulting in six random differential equations.
We rewrite our general neural mass formulation as The neural masses in a single ZJR cortical column have different properties and are characterized by the parameter matrix Expressing the well-known ZJR cortical column with our elaborate formulation might seem a trivial exercise.However, it does prove that model includes the simpler cases.When dealing with more complex models, such as those described in the next section, the advantages of the general model might become more evident.

Population-Level model as the result of the interactions between several ZJR-cortical columns
Some literature explores a double-column NMM (57,70,71).We thought it instructive to test our NMM with a larger number of cortical columns a model we call "Population-Level".We now need to label each ZJR cortical column with a superscript ( u U  ) and a subscript to any quantity that refer to the three types of neural mass in each column (i.e.,

 
Pyr, Inh,Ste j  ).For example, the system parameter matrix for the unit u is ( ) ( ) ( ) The complete system parameters are contained in the matrix ( ) ( ) ( ) , , This time, the probability density of transmission delays is ( )

Ste
).The red arrows describe the intrinsic connections between neural masses, while the black arrows represent the inputs to the cortical column affecting the Pyramidal neural mass and the outputs via the Stellate neural mass.
where the submatrices s C and l C define the short-range and the long-range intercolumn connectivity, respectively.
The scalar floor 6 Nu Ns represents the number of sample points for the short-range connectivity and Nu is the total number of cortical columns involved in the model with , i j U  .
To model short-range connections we use an exponential function (31,59), defined as where b is a constant.We determine the long-range connections randomly, which could be extracted from diffusion analysis.
Here, we illustrate with an example of five ZJR cortical columns (i.e.,

Definition of Connectivity Matrix topologies 0 K and conduction delays for NMM simulations
We explore the effect on neural dynamics of three different types of Connectivity Matrix topologies 0 K and delays.
For this purpose, we assume that the neural masses are distributed on a ring in two dimensions.
• For the first Connectivity Matrix topology, we specified short-range connectivity between nearest-neighbor nodes on the ring.Large-range connections were not considered, i.e., 0 l = C (See Figures 5A1, A2).We call it NN network.
• We keep short-range connections to neighbors in the second case, but we define sparse large-range connectivities (See Figures 5B1, B2).We call it SW network.
We call it FC network.For all simulations, delays within a column are instantaneous.We explored the three different types of delay tensors D described in section 2.1 and illustrated in Figure 6:

Simulation of a single ZJR cortical column
The classical ZJR NMM of a single cortical column comprises three neural masses (See Figure 7) with a delay tensor 1 D corresponding to zero delays.The simulation used a delta of 1 msec, and a burn-in period of 0.5 sec.After the burnout, a segment of 1 sec was retained for the analyses.In what follows, we analyze the output the Pyr neural mass.
At this level, we tested 3 different scenarios for our mesoscale simulations.Each simulation scenario used the values for the dynamic noise, as shown in Table 5.The first (determinstic) scenario has a relatively high mean and rather small standard deviation, i.e.

30, 1
 == (Figures 7B1, B2, and B3), and the third scenario retains the small standard deviation but also has a small value for the mean, i.e.

3, 1
 == (Figures 7C1, C2 and C3). Figure 7 displays an example of the 3 simulation scenarios.Our analyses are carried out to simulate neural mass output Pyr time series and estimate the corresponding spectrum by the multi-taper method from the Chronux Matlab package (73,74).The rightmost column shows the 3D of phase plot for the output of the 3 neural masses: Pyr, Inh, Ste.
The simulated Pyr time series for all three simulation scenarios are smoothly periodic.As we can be seen, the first two simulation scenarios are pretty showing only light perturbations in the spectrum and into the neural dynamics of the second simulation (See Figures 7B2 and B3).These perturbations are due to introducing a small noise level in the NMM.
We observe that for the first two scenarios (i.e.,

30, 0
 == and 30, 1  == ), the spectrum loses the harmonics for frequencies higher than 60Hz , while the spectrum preserves more harmonics for the third scenario (i.e.

3, 1
 == ).As was expected, the dynamics of the three simulation scenarios all converge to a limit cycle in the state space plot, regardless of the dynamical noise values explored.

Simulation of several ZJR cortical columns
We then tested our toolbox with an extended neural mass network (eZJR NMM).The number of cortical columns was increased to 1000 columns = 3000 neural masses.We restricted attention hereafter to the first (deterministic) simulation scenario 30, 0  == .Here we explored the effect of changing two aspects of the Connectome Tensor.
1. We investigate whether the eZJR NMM is affected by the topology of the connectivity matrix topology ( 0 K ).
2. We examine the influence of different types of delay tensor distributions (D ) on the dynamics of the ZJR NMM.

K ) on the eZJR NMM
For these simulations, we fix the delay tensor 3 D (Figure 6C).We consider three different topologies of connections.
• The NN network: This type of network comprises only short-range connectivity between nearest-neighbor nodes without large-range connections (Figures 5A1, A2) • The SW network: A "small world" type of networking retains short-range connections to neighbors and includes sparse large-range connectivities (Figures 5B1, B2) • The FC network: A fully connected graph (Figures 5C1, C2) The analysis of the output eZJR model is similar to the one made for a single ZJR cortical column.However, instead of looking at the output of a single ZJR NMM, we study the average behavior of the NMM ensemble.Therefore we generate a sample 3000 time series, one for each NMM.The results are in Figure 8.Here, the leftmost column is the average time series of 1000 Pyr NMM (Figures 8A1, B1, and C1).The middle panel (Figures 8A2, B2, and C2) shows the corresponding power spectrum of the average Pyr NMM outputs.Finally, the rightmost column is the phase plot of the eZJR NMM (Figures 8A3, B3, and C3), each axis reflecting the average of the respective types of neural masses: Pyr, Inh, and Ste.
We got periodic time series and similar spectra, showing all the harmonics independently of the Connectivity Matrix topology (Figures 8A1, A2, B1, B2, C1, and C2).In the phase portrait of our eZJR NMM, the dynamics correspond to stable closed curves but are a bit different among them (Figures 8A3, B3, and C3).As we can see, the FC network dynamics is the simplest and the most similar to a limit cycle(Figure 8C3), while the NN network is the most complex with an extra protuberance (Figure 8A3).The spectra show quite differences among them, which could be associated with the linear spectra used to make the approximation.

Influence of different types of delay tensors (D ) in the eZJR NMM dynamics
For the second study at the Population Level, we focus on the influence of time delays into the deterministic ZJR NMM simulation.We introduce three different distributions to represent the activity propagation modeled by a delay tensor D .We investigate the influence of the three delay tensors: • The instantaneous delay tensor ( 1 D ), • The synchronized delay tensor ( 2 D ), • The distributed delay tensor ( 3

D ).
We defined them in Section 2.1 and after mentioning them in Section 3. To carry out these simulations, we focused the SW network.
As we illustrate in Figure 10, the simulated time series still behaves smoothly periodically; perhaps we model a considerable number of unsynchronized ZJR cortical columns (Figure 9A1, B1, and C1).Their spectrums are relatively similar with minor changes, preserving the alpha peak and the harmonics (Figure 9A2, B2, and C2).The trajectories in the phase portrait, although different, are relatively similar, mainly for the first two delay tensors (

Discussion
As mentioned before, despite neural mass models (NMM) being a mainstay of current neuroscience modeling and analysis, current toolboxes suffer from limitations due to the choices of model formulations (9,10,(12)(13)(14)(15).The most popular methods are based on tightly coupled systems of simultaneous differential equations whose integration quickly becomes computationally prohibitive as each neural mass's complexity or the system's size to simulate increases (17)(18)(19)(20)(21). Additionally, usual differential equation techniques that seem popular have two fundamental limitations.i) They use non-optimal methods for integrating stochastic or random systems of differential equations.Some of these methods change the original continuous-time system's dynamical properties (attractors).ii) There is no optimal way to deal with transmission delays between neural masses.There isn't any solution to incorporate into neural mass modeling distributed delays, which is the natural formulation considering the diversity of axonal myelination and diameters.
We leverage a tensor formulation of NMM to solve all these issues:  D ) and Bottom: we describe the model using a distributed delay tensor ( 3 D ).We make an analysis for each type of delay tensor, based on the generation of the simulated time series (A1, B1 and C1), the spectrum (A2, B2 and C2) and the oscillatory dynamics (A3, B3 and C3).
• Recognizing that neural transmission occurs at a finite speed, we solve the system of equations with a sufficiently small integration step to ensure that input to each neural mass from other masses can be calculated before solving the neural mass equation.This approach uncouples the integration of each neural mass, with significant simplifications in the solution of the systems of differential equations (which are now low dimensional) (20,25,26).• In turn, the integration of each neural mass may now be optimized.We employ the local linearization method, which preserves the dynamical properties of nonlinear systems (27).By decoupling the solution of the equations for each neural mass, we overcome one of the prime difficulties of the LL technique, which involves using the exponential matrix of the Jacobian of the system.Without decupling, this operation becomes computationally demanding at an explosive rate with an increasing number of neural masses.Decoupling the solutions of NMM makes this integration feasible.Simple NMM models such as the Zetterberg-Jansen and Rit (ZJR) (45,57) or Wilson and Cowan (75,76) models can be solved symbolically (before the simulation) to great computational advantage.• The final ingredient to our formulation accrues from stating the model, not in terms of firing rates but rather in terms of postsynaptic potentials (31).The sigmoid function is applied to the NMM input to produce its output.The input is a linear combination of the past outputs of all other neural masses for a finite time window.It follows readily that the input of the neural mass can be stated as an operation between the Connectome Tensor of delays and the previous output of all other neural masses.• The Connectome Tensor efficiently encodes any realistic delays between neural masses, including distributed ones, without recourse to ad hoc solutions such as prolonging the postsynaptic potentials.
For single neural ZJR neural mass, the results produced by our toolbox are identical to previous results (45,48,57,77,78) with the same system parameters.The simulation with 1000 cortical columns would have been infeasible with other approaches.The simulations presented in this paper are relatively simple.However, it is striking that different network topologies and delay distributions produce roughly similar results, though essential nuances are evident in the results of the simulations as we tweak the NMM parameters.We plan to use this toolbox to carry out more extensive and realistic modeling, incorporating more biophysical and physiological details.
Two areas need to be studied further.Tensor computations have been optimized for high-performance computer systems.We must modify our code to fully take advantage of the tensor formulation to speed up the computations.This improvement is especially needed when inserting our algorithm into the solution of inverse problems in which simulations must be repeated many times to estimate system parameters.

Code Availability
The procedures described here are in a toolbox that might enable more realistic large-scale brain modeling, allowing a high-performance simulation of neural dynamic behavior via fast and efficient implementation of the essential core of calculations.Our toolbox and codes are publicly available via Github (https://github.com/CCCmembers/Neural_Modelling).
To calculate the time series spectra, we use the Chronux package, an open-source software package for analyzing neural data, available via http://chronux.org/.

Author Contributions
The contribution of AGM to the presented results was the design of the research, the development of the theory, mathematical demonstration, design of methods, paper preparation and implementation, and the development of the toolbox with DPL and AAG.ML and YW collected the essential references in this field.PVS and MBV senior authors supervised the process and did the final revision.PVS introduced the theoretical background, helped write the article, and guided this work.All authors contributed to the article and approved the submitted version.

Figure 1 :
Figure 1: Connectome tensor of our Neural Mass Model.A) shows the physiological relation between two group of neural masses, where axons, emitter and receiver neuons are involved.The base of our NMM is a tensorial formulation of this phenomen, which we describe in B).Figure created with BioRender.com.

d, 1 D and 2 D ) 0 K
A time delay between an emitter and a receiver node e r Length of the fiber system connecting an emitter and a receiver node  continuous or discrete time (Examples explored in this work 0 D Connectivity matrix that does not explicitly model neural transmission delays ,

Table 6 : 2 . 1
Set of parameters to configure the Zetterberg-Jansen-Rit (ZJR) NMM General model formulation Our general neural mass model is expressed in the state-space form: State equation:

.
as the set of possible transmission delays, with max  the maximum delay of information transfer between two neural masses.We also introduce the delay tensor ( ) We can now define the continuous-time Connectome Tensor ( same dimensions as D via the Hadamard product:

a
Dirac-delta centered at 0. Then = 0 K K and max 0  = , simplifying to the usual no-delay connectivity model.

.
We now discretize the Connectome Tensor as: is the probability density of transmission delays for each delay model 0,1, 2 i = defined above, and

jh
is the maximum amplitude of the PSP and j  is the time constant of the membrane involved in the model, which depends on the rate constants of passive membrane and other spatially distributed delays in the dendritic tree(58).
refers to the nonlinear sigmoid (the next potential-to rate) function that determines the pulse density and transforms the average membrane potential of a neural population into an average firing rate(57,65).It is defined as

Figure 2 :
Figure 2: This figure describes the dynamics in a single neural mass.The random inputs to the neural mass are collected in the term ( ) ( ) j j j j t w t  =   +   see Figure BB, which contains the model parameters of the three neural masses.The system inputs act on the Stellate cells, and the neural activity comes out from the Pyramidal cells.See Table4for further description.For a single ZJR cortical column, the delay tensor 0 D is modeled by ( is, only with local intra-columnar connections, and we assume that the transmission time is faster than the sampling rate.Therefore, the Connectome Tensor reduces to weights of the directed connections between the neural masses, which are well established in the literature to generate EEG-like activity(57,69).Furthermore, the synaptic connectivity between Excitatory (positive weights) and Inhibitory (negative weights) cells is shown below

,
which means that the delay tensor corresponds to 2 D .We build the Connectivity matrix 0 K manually, satisfying that

Figure 3 :
Figure 3: This a scheme of a single ZJR cortical column that include three neural masses (Pyramidal: Pyr , Inhibitory: Inh and Stellate:

5
Nu = ) how the nodes are connected Therefore, once the connectivity matrix is established, we can compute the CT as is the time lags for two connected neural masses.Afterward, we calculated the output of the model at the Population Level as modeling of a single cortical column is a particular case of the Population Level, where 0 K coincides with the intrinsic connections in a cortical column, i.e.,

Figure 4 :
Figure 4: This figure shows how five ZJR cortical columns are connected, where the double arrows in red and black color symbolize the shortrange and the long-range connections, respectively.

• 3 D
(distributed delay tensor), the gamma distribution function for axonal delays proposed by Nunez(51,72), where the probability density function for corticocortical propagation peaks is about 6 distribution estimated to be about 3 m / s (Figure6C).

Figure 6 : 1 D 2 D 5 . 3 D
Figure 6: We explore three types of propagation velocity to determine the conduction delays and here we show how is the function behavior.The first way (A) is using 1 D , for all connections a Dirac delta function centered in 0 , while the second one (B)

Figure 5 :
Figure 5: An example of the Connectivity Matrix topology for 15 ZJR cortical columns placed on a ring.The first row displays a circular graph visualization (Created with circularGraph), while the second one illustrates the sparsity pattern of the Connectivity Matrix.The first column (A1, A2) corresponds to neaerest-neibhbor connections, without large-range connections.The second column (B1, B2) Adds sparse large-range connections to the nearest-neighbor ones.The third column (C1, C2) is a full connected topology.

Figure 7 :
Figure 7: This is a three-panel figure of a single ZJR cortical column simulation for three examples of dynamic noise values.Results are shown in the left and middle columns for the model output (Pyr).The right column shows the state spce represnatation of the activity of the 3 neural masses (Pyr, Inh, Ste).Top row: simuaton without noise, i.e., 30, 0   = = , Middle row : the simulation with a small increase in the standard

1 D and 2 D
) employed.They are stable closed curves showing up the stability of the deterministic system.

Figure 8 :
Figure 8: Three-panel figure for the connectivity analysis under three types of connections throughout 1000 ZJR cortical columns.Top: correspond to NN network, Middle: represent the SW network and Bottom: describes the FC topology.We make a study for each type of connectivity, based on the generation of the simulated time series (A1, B1 and C1), the spectrum (A2, B2, C2) and the oscillatory dynamics (A3, B3, C3).

Table 1 :
Structure Types

Table 4 :
Set of variables used in our Neural Mass Model
Ste tx ,