Neuroinformatic tool to study high dimensional dynamics with distributed delays in Neural Mass Models
A. González-Mitjans, D. Paz-Linares, A. Areces-Gonzalez, M. Li, Y. Wang, ML. Bringas-Vega, P.A Valdés-Sosa
CComputational tool to study high dimensional dynamic in NMM
González-Mitjans A. , Paz-Linares D. , Areces-Gonzalez A. , Li M. , Wang Y. , Bringas-Vega M.L. and Valdés-Sosa P.A.
1. The Clinical Hospital of Chengdu Brain Science Institute, MOE Key Lab for Neuroinformation, University of Electronic Science and Technology of China, Sichuan, Chengdu, China 2. Department of Mathematics, University of Havana, Havana, Cuba 3. Department of Neuroinformatic, Cuban Neuroscience Center, Havana, Cuba 4. Department of Informatics, University of Pinar del Rio, Pinar del Rio, Cuba
Abstract:
Neuroscience has shown great progress in recent years. Several of the theoretical bases have arisen from the examination of dynamic systems, using Neural Mass Models (NMMs). Due to the largescale brain dynamics of NMMs and the difficulty of studying nonlinear systems, the local linearization approach to discretize the state equation was used via an algebraic formulation, as it intervenes favorably in the speed and efficiency of numerical integration. To study the space-time organization of the brain and generate more complex dynamics, three structural levels (cortical unit, population and system) were defined and assumed, in which the new assumed representation for conduction delays and new ways of connecting were defined. This is a new time-delay NMM, which can simulate several types of EEG activities since kinetics information was considered at three levels of complexity. Results obtained in this analysis provide additional theoretical foundations and indicate specific characteristics for understanding neurodynamic.
Keywords:
Neural Mass Model, time-delay, connectivity matrix, brain dynamic, kinetics information, local linearization method.
Introduction
Representing real-life practical problems through mathematical models is a challenging task. One of the factors to consider in neuronal modeling is the terms of experimental validation and the use of physical units to represent them. Because of this, many mathematical models have been developed to represent and describe the human brain, regardless of the level of abstraction and the assumptions that have to be made for its formulation. A mathematical description of a neuron that predicts its dynamics and characterizes its properties, is known as a Neural Mass Model (NMM), which are represented by a set of nonlinear differential equations and often have stochastic information. One of the oldest and most recognized models is the Hodgkin-Huxley model (1963), which is defined by a system of nonlinear ordinary differential equations that reflect the propagation and emergence of action potentials. Electric potentials are a consequence of spiking neurons. The nerve cells that present this dynamic type are of greater importance in the biological study of neuron models. There is a large amount of literature on the implementation and discussion of all these models. Due to the need for more realistic representations of the brain as a whole, these models have a high dimensionality, which further complicate its treatment. Certain neuroscientists have considered this aspect and there are several software and packages such as GENESIS and NEURON, which are used to silico modeling of realistic neurons [1] and others as The Virtual Brain (TVB), The Big Brain (TBB) and the Blue Brain (BB) project, which aim to construct a biophysically detailed simulation of a cortical column on supercomputer with a high anatomical definition. The main point in these tools is that detailed neuron descriptions require methods computationally expensive and without a fine time scale. Therefore, fast developments of a high-performance computing system are necessary. Based on this, the research content has covered the implementation of a brand-new procedure [2] , which instead of using Stochastic Differential Equations (SDE) relies on Random Delay Differential Equation (RDDE). This enables the solution of the integration of the neural masses step by step for more efficiency and at a really fine time scale. It would appear if one uses a really fine times scale (impossible in the other formulation) then what would happen is would be a very long time. On the contrary, if one uses a very fine time step, the events of the neuronal simulations become unconnected, since then inputs for any given neural mass are delayed and are fixed values the generation of the postsynaptic potentials and other states. The solution of this set of RDDE become uncoupled for each neural mass and they can be solved explicitly with a computer algebraical program as demonstrated in the paper. In addition to the numerical implementation of the ingenious contribution created by Valdes-Sosa, a simulation of the human brain was carried out, including three dynamic levels. A brand-new analysis for the conduction delay between cortical units at the population and system-level was thought of, wherever the distribution used by certain researchers to model the time-delay in the system is modified. A rapid simulation was performed with a high number of ortical columns, taking into account the kinetics information of the brain, which allows a representation of all frequency bands.
Model Description
Jansen and Rit Neural Mass Model
NMMs describe by mathematical formulation neural activity at the mesoscopic level. One of these most famous models is the one formulated by Jansen and Rit. The essence of this model is modeling through three coupled nonlinear differential equations of the second order, which are transformed into a first-order system with six differential nonlinear equations. The Jansen and Rit model performs a transformation of the action potential into the average of the postsynaptic potential, creating a delay in the synaptic transmission, which is essential to generate the oscillatory activity [3] – [6] . Many have been the studies made on Jansen and Rit and several changes had been done in its formulation. Single Cortical Unit
This time, we deal with a new interpretation of the Jansen and Rit model known from Valdes, 2009. This NMM contains three neural masses in each cortical unit. The pyramidal (Pyr) cells, inhibitory (Inh) and stellate (Ste) layers are used to generate de average postsynaptic membrane potentials. Valdes, in his study to the classical Jansen and Rit NMM, modified the classical model adding an algebraic term to make calculations easier and faster. The current study maintains this ingenious change and provides a phase state analysis to study the behavior of two or more connected Jansen and Rit neural populations. Some research has been done about it taking a fixed value as the input and considering a time-delay as a Dirac- distribution. Here, a new interpretation of the delays in the model is examined and the kinetics information in the brain regions. A single unit is given, in this model, by the interaction between three neurons (or layers), two excitatory and one inhibitory. This formulation is considered in this work as the most basic level of brain dynamics. The state-space of the model and its internal dynamic between layers into a cortical unit are shown in Figure 1. The pyramidal neurons in this model receive feedback from the stellate and inhibitory neurons. The input to a neural mass is transformed from an average pulse density to an average potential membrane, via postsynaptic potential function. Then, the output of the system will be a convolution of the input and the impulse response function. The impulse response function of the model (or the first-order kernel) is interpreted as the characteristic PSP elicited by a single incoming spike given by a Dirac- function as follows
0( ) 0 0 l tl ll H te th t t − = (1) where l H is the maximum amplitude of the PSP and l lump the characteristic decays transmission together, i.e. the time constant of the membrane with the different layers
1, 2, 3 l = involved in the model. l is a lumped representation of the sum of the rate constants of passive membrane and other spatially distributed delays in the dendritic tree and satisfies that l l = [7] . The values of these parameters depend if the EPSP or the IPSP [8] – [10] . The maximum value of the EPSP is reached faster than the maximum of the IPSP, but the latter has a greater potential. There is also a static nonlinear sigmoid (the next potential-to rate) function, which transforms the average membrane potential of a neural population into an average firing rate [11] – [13]
0( )
2( ) 1 r v v eS v e − = + (2) where e is the maximum average firing rate of the neural populations, v is the PSP corresponding to e and r is the slop of the sigmoid . The current NMM is formulated trough an algebraic equation [2] , which allows the use of more efficient
Figure 1: Jansen and Rit Neural Mass Model scheme for a single cortical unit. At this level the interaction between excitatory and inhibitory layers is studied. ools. The model is formulated as a set of first-order nonlinear random differential equations with the state vector ll l yx = x . The connections and the information about conduction delays was considered as follows ( ) ( )( ) [ ( ) ( ( ))]2 ( ) ( ) l ll l l l l l ll l l l y t x tx t H w t S z tx t y t = = + + − − (3) ( ) ( ) l l z t y t = K (4) The main variables of the model are represented by l y which is the output of the PSP and also their derivatives were considered l x for each layer
1, 2, 3 l = . Where the relation between the values of the layer l and the name of the neural mass is showed below lll = = = The random properties of the system are driven by a noise process ( ) l l l w t + , which is renamed by l p , with
1, 2, 3 l = as is showed in Figure 1. It represents the first derivative of a Brownian noise process with drift and infinitesimal variance l where ( ) l w t is the first derivative of a unitary Brownian process. The Algebraic Differential Equation (ADE) (4) involves the variable ( ) z t , which is a linear combination of the PSP and the values of the synaptic connections between the neural mass
1, 2, 3 l = . Moreover, the output MEG/EEG signals are represented by the first component of this term (i.e. ( ) z t ). K defines a connectivity tensor with the conecctions and time-delay information. It is a combination of the classic connectivity matrix l l c = C , which contains the weights of the directed connections between two layers , , , 1, 2, 3 l l l l l l → = , and an exponential distrubution to represents the time dealys in the model. More details about this term will be exposed in next sections. Table 1: Standard values of Jansen and Rit NMM parameters
Parameter Physiological interpretation Value l H Maximum amplitude of the PSP for excitatory and inhibitory neural masses respectively l Intrinsic lump parameters for excitatory and inhibitory neural masses respectively
10 ms, 50 ms l Lump the characteristic decays transmission together for excitatory and inhibitory neural masses respectively l l = e Maximum average firing rate of the neural populations v The value of the potential r Slope of the sigmoid at v -1 c Synaptic connections c , c Average number of synaptic contacts in the excitatory feedback loop
13 31 c c = = , c , c Average number of synaptic contacts in the slow feedback inhibitory loop
12 21 c c = − = If l l = , then l l c = as is illustrated in the formulation (5). Its off-diagonal values are expressed as fractions of a constant c = , based on empirical knowledge of the EEG alpha-like activity and estimation made from references [12][14] .
12 1321 31 c c c cc c c c = = = = (5) The connectivity matrix used in the model is generated according to the connections shown in Figure 1. n the other hand, if it refers to a hyperpolarization then the ions moving across the membrane are negatively charged and correspond to an Inhibitory Postsynaptic Potential (IPSP) [15][16] . This theory was perfectly exposed by Hodgkin and Huxley [17] . Furthermore, the synaptic connectivity between neuron masses is positive except when it is from pyramidal to inhibitory cell (i.e. l l = → = ), in this case, it is negative. − = C (6) Model parameter values variables and functions of Jansen and Rit single column are given in Table 1 [12][13] . The brain is much more complex than we can imagine, so a more detailed interpretation is necessary. In the next section, it scales a little more , taking the examination at the population and system level, which are defined in this work. Model extension to complex levels: Population and System Levels
In this study, the interaction between two or several cortical units with similar characteristics was called population level. A schematic of the interaction between two cortical units is shown in Figure 2. At this level, the response function is like that of a single cortical column, but the but the number of cortical units involved in the system is also considered
0( ) 0 0 lu tlu lulu H te tR t t − = (7) where lu H and lu are the same parameter described in the previous section with
1, 2, 3 l = are the neural masses and
1, , u Nu = represents the cortical unit and Nu is the total number of cortical units involved in the model. Moreover, the sigmoid function is exactly the same and the NMM includes information about the cortical unit as is described in the equations (8) and (9). ( ) ( )( ) [ ( ) ( ( ))]2 ( ) ( ) lu lulu lu lu lu lu lu lulu lu lu lu y t x tx t H w t S z tx t y t = = + + − − (8) ( ) ( ) lu lu z t y t = K (9) Here, the variables and parameters involved were exposed before for each
1, 2, 3 l = and now, the cortical units
1, , u Nu = are also considered. The connectivity between cortical units is given by the tensor K with the time-delay information and coefficients of the connections between cortical units , , 1, , u u u u Nu = . This connectivity tensor is computes as = K D K (10) K is a matrix with the inter cortical unit connections and the connection strength tensor D contains the weight-delay values. Figure 2: Interaction between cortical units defined the Population level o build the system level, information from several neural populations is needed. To the information obtained via population level is added neural kinetics features, those allow the simulation of the relationship between different brain regions with different rhythms [7] . Taking this in account similar procedures were made to simulate population and system levels. Since the NMM is formulated by a stochastic nonlinear system, the analysis of this model is more complicated. That's why a linearization of the model was performed, using a method developed by Ozaki and colleagues.
Local linearization method
The convection of the traditional Jansen and Rit model to a system of algebraic differential equations was an important transformation to counteract out the problems of scale integration and facilitate the analysis of very large networks of neural masses with an efficient resolution [2] . Here, the LLM of RDEs on the model was employed [18] . This method for the single cortical unit will be exposed below. As is known, the manipulation of nonlinear systems is more complex, so it is widely important to find a linear equivalence. el to a system of algebraic differential equations was an important transformation to counteract out the problems of scale integration and facilitate the analysis of very large networks of neural masses with an efficient resolution. Here, the LLM of RDEs on the model was employed. This method for the single cortical unit will be exposed below. As is known, the manipulation of nonlinear systems is more complex, so it is widely important to find a linear equivalence. Consider the state space of the model for a single cortical unit, given by equations (3) and (4). Let ( ) l p t , ( ) l t x be the stochastic term and the state vector of the model respectively, where
1, 2, 3 l = for each , t t T and f characterizes the k -dimensional RDE of the NMM as follows ( ) ( ) ( ) ( ) ( ) , l l ll l d t f t p t dtt == x xx x (11) As this is a system of nonlinear stochastic differential equations, it can be solved numerically via LLM with a given step size for integration, i.e. n t t nh = + and n = . Then, f J hn nl l e + = + x x r L (12) where Tk + r = , x2 , 0 k k = L I , with k = and f J is defined by blocks as is shown in equation (1) in the Appendix section, where l x f J and l p f J indicate Jacobian operators of the function f concerning each of its components. Algebraic derivation was performed by the MATLAB Symbolic Toolbox; a detailed procedure is given on a live script. The results of these operations are shown in equation (13), where , l l H are the PSP parameters described in the previous section with
1, 2, 3 l = , ( ) ( ) l dS z t and is the derivative of the sigmoid function. ( ) ( ) ll l lp l l l l l H dS z t H = − − = x ff JJ (13) Furthermore, the full expression of f J is a ( ) matrix, obtained after substitute the Jacobian matrices of f (See equation (2) in the Appendix section). In the next step of the LLM, after to apply the exponential matrix, a regrouping of terms was made, which simplifies and facilitates the compression and manipulation of the following numerical expression ( ) , , , l t h l t l t + = + + x x I A BE (16) where A and B are matrices with the PSP constants for each
1, 2, 3 l = neural mass, I is the identity matrix and , l t E contains the random inputs to the system and the sigmoidal function (See equation (3) in the Appendix section) [2] . The terms diff ( ) l p t and diff ( ) l z t in the , l t E matrix, represent the difference between the values of ( ) l p t and ( ) l z t in the current and the next step respectively. To make the work easier, an auxiliary three-dimensional matrix of dimension ( ) was created to carry out the algebraic numerical integration, which contains the PSP integration factors for each
1, 2, 3 l = . This new term grouping has an advantage, which is that it allows you to create a multidimensional array (also interpreted as a tensor) that will not change. The matrices A and B respectively for each layer
1, 2, 3 l = are constants egardless of the level being analyzed, which favors numerical integration. At system level, their dimension increase, containing the initial matrices A and B for each population (gamma, beta, alpha, theta and delta). While the arrays containing the PSP info do not change, the way the MEG/EEG output is generated changes depending on the level being analyzed, which influences the value of , l t E or , u t E according to the level studied. But it is necessary first to define how the connectivity tensor will be created and how the conduction delays will be assumed, compute ( ) l z t (or ( ) u z t ) and then , l t E (or , u t E ), without loss of generality. Testing the influence of different time-delays on the oscillatory behavior of the model
We examined three types of delays, two of them study previously. First, we test with a fixed value, which represent the “first generation of NMM s ” via a Dirac- function centered at the origin and the “second generation of NMMs ” given by a translated Dirac- . Some researchers have focused on finding a good approximate to represent the time delay suffered by the neural transmission of information [19] – [21] . Despite including delays, these models still had problems. The conduction-delays considered until now make the temporal dimension not well represented. Time lag has been explored by previous researchers for a Dirac- centered at the origin, or translated. [4][22]– [25] . It can be written as a sum of two points connections, given for two different layers , , 1, 2,3 l l l l = as follows ( ) ( ) l l l l l ll z t c y t = = − (17) where l l represents the conduction delay estimated as a multiplication of the average path length between voxels (nodes), and conduction velocity is
10 m s [2][7][12][13] . As many of the natural phenomena have an exponential evolution in time, we studied how the model behaves with a distribution of this type. We found the model is sensitive to changes in time delays, which is evidenced in the stability of the system. Therefore, the way in which the delays are represented is very important. As evidenced in Figure 3 the new interpretation for driving delays implemented in this work provides greater stable oscillations for the NMM.
New Interpretation for Weight-Delays and Connections
To determine how the output of the neural model will be estimated, first it is necessary to propose how to compute the procedures at these levels and define the new interpretation for the conduction delays. These assumptions will provide greater similarity to real neural activity. Consider a tensor u u D , with weight-delays between neurons and was defined as an exponential Figure 3: Brain dynamic changes according to time-delays istribution. The matrix u u K contains connection information between cortical units involved , , 1, , u u u u Nu = . Therefore, the connectivity tensor for the population and system level with delays information and coefficients of the connections between the cortical unit u and u was defined as u u u u u u = D Κ K (19) Figure 4: Tensor representation of connections and weight delays in NMMs
In a similar way, a connectivity tensor was calculated for the simplest level (cortical unit level), but the dependency relationships were established between neurons (Pyr, Inh and Ste). The procedure to compute this multidimensional matrix was expressed by a tensorial appreciation. Then, the output of the model at the population and system levels is giving by ( ) u t z , with , , 1, , u u u u Nu = . It may be calculated through the double sum of the connectivity tensor and the EPSP emitted from other units added to the MEG/EEG output l z of the cortical unit level ( ) Nu Nu l u u u u uu z z y t = = = + − Κ (20) where
1, , N refers to the time lags between two connected cortical units , , 1, , u u u u Nu = . The difference between population and system outcomes is that to obtain the system output, the kinetics of each population (from delta to gamma rhythms) were considered. Taking this into account the matrix , u t ε is given by ( ) ( ) ( ) ( ) ( ) ( ) , diff diff ( ) u uu t u u u p t S z tp t z t dS z t + = + E (21) here, diff ( ) u p t and diff ( ) u z t are the difference between the current and the next step of the values ( ) u p t and ( ) u z t respectively, with
1, , u Nu = . The random term ( ) u p t in this case constitutes an extension concerning the previous one in terms of dimensionalities. Time-delays was simulated as an exponential distribution between two different units and populations. To generate the center of the distribution, it was necessary to first select a frequency (in Hz ) that would be present in the cortical unit and characterize the frequency of the population. The values of a lumped representation of the rate constants of the passive membrane and other spatially distributed delays in the dendritic tree ( , 1, 2, 3 l l = ) to generate the different frequencies bands were estimated from reference [7] . Therefore, five populations will be created, each characterized by a different frequency band (delta, theta, alpha, beta, and gamma) and the neural dynamics between them will define the system level. The excitatory layers are given by the Pyr ( l = ) and the Ste ( l = ) layers and the inhibitory, as its name implies by the Inh ( l = ) layer. Then, it will be considered , e i for the excitatory and inhibitory layers, respectively. These parameters were obtained through a comparative and analogous estimation, taking data form previous works, and extrapolating it to provide the estimation for the current research. The estimated values of lump parameters used to obtain oscillatory behavior are shown in Table 2. Table 2: Estimation of the excitatory and inhibitory intrinsic lump parameters
Frequency bands e i gamma (>30 Hz) 1 ms 4.8 ms beta (12-30 Hz) 4.5 ms 10.7 ms alpha (8-12 Hz) 10 ms 20 ms theta (4-8 Hz) 15 ms 40.7 ms delta (1-4 Hz) 20 ms 52 ms There is a constant ratio for the values of both, excitatory and inhibitory layers, which was very useful to estimate the intrinsic lump parameters and it is formulated as ll lH l == = (22) For the simulation at the unit and population level, it is only necessary to select one of the frequencies resented in Table 2, according to the interests of the researcher. While generating the dynamics at the system level it will be necessary to consider all the frequencies bands, showing the differences between populations with different frequencies and will make the simulation more realistic. To compute a tensor with the weight-delay information D , were considered the center of the delay distribution and the width of the delay distribution , then exp t − = − D (23) These parameters were obtained through the process to create synthetic directed connectivity and delays, which internally needs the frequency that will be used in the cortical unit and population simulations or a set of frequencies for the system level. As for the connectivity of the model, at the cortical unit level was used the classic connectivity matrix estimated by Jansen and Rit, which was exposed in the formulation (6). The connections at the population and system levels in the NMM was implemented by a new formulation. First was estimated a matrix K using spectral factorization. Consider the following expression ( ) ( ) † 3 = − − Q I K I K (24) where the sigma (variance) of the Wiener process for Stellate layer, which constitutes the input layer of the model and Q is the synthetic directed connectivity reconstruction generated via the Hermitian Gaussian Graphical Model to construct a multivariate complex Gaussian and sampling them [26] . Let be auxiliary matrices calculated from singular value decomposition and satisfy that ( ) ( )( ) ( ) †3 33 == = † † Q UDUQ U D U DQ U D U D (25) Therefore, ( ) = - K I U D implies that ( ) ( ) t = − F K K and ( ) t − is a Kronecker function. represents the intrinsic frequency response used in the simulation. From this procedure the values of the parameters , , and K are obtained, which are the center of the delay distribution, the maximum of the lag distribution, the width of the delay distribution and the information about the connection between the cortical units (or populations) involved. For the population analysis, the connections are defined from , , u u u u → → with , , 1, , u u u u Nu = and for the system level, the connections are given two to two between the five populations generated by different frequencies. These results are used to compute the connectivity tensor Κ which is widely important for generating the input to the system ( ) u z t , with , 1, , u u u Nu = . Results and Discussion
In the current work was simulated electroencephalography (EEG) data generated by a new formulation of the Jansen and Rit Neural Mass Model (NMM). Due to the existence of the time delay needed to propagate the actual signals among neurons, a new time-delay NMM and tools to study its dynamical behavior were developed rather than maintaining the tendency of previous researchers to employ a Dirac distribution. It was established as an exponential distribution, effectively simulating the conduction delays of the signals across the axon and its advantage was proven. Conduction delays and connections were combined on a clearly tensor implementation, making the temporal dimension well represented. The study of NMM in a single cortical unit was extended, developing a brain simulator at three levels, including kinetics information (gamma, beta, alpha, theta, delta) sown in Figure 5 and Figure 6. Through the simulated EEG recordings for each rhythm using the new form assumed to consider conduction delays, it was found that there is a greater similarity to the signals emitted during real EEG analysis and they are consistent with the values of these frequencies and their neural dynamics. The brain simulations reach a limit cycle behavior depend on decreasing the frequency. For larger cortical units (1000 units) stable simulations were obtained that could produce many types of EEG frequency spectra and stable behavior was found also mixing units with different characteristic frequencies (See Figure 5). hese findings can indicate specific characteristics for each frequency bands, which coincide with the idea developed by Lopez da Silva about the importance of the oscillations emitted by lower frequencies to establish a functional bias in populations with a large number of neurons. It can also be of great relevance in obtaining specific response information in resting state and on neuronal sets interacting under a particular task. The implementation of the computational platform developed to carry out these studies were fast and efficient despite the high dimensionality of the problem and the inclusion of kinetics information. This was made possible by an algebraic reformulation of the classic Jansen and Rit model, which speeds up the speed of the program and makes calculations for numerical integration easier. This is a very important aspect as many of the researchers in this field are using traditional methods to integrate numerically, such as Runge Kutta and Euler, which have a high computational cost. The algebraic numerical integration was done in Live Script integrated into MATLAB as a tool to be disseminated.
Figure 5: Brain simulation at cortical unit and population levels. Each row represents different frequency bands for the two simulated levels. The first column of each level contains an EEG simulated recording and a time-frequency map, while the second one describes the oscillatory behavior of neurons and units (1000) at unit and population levels, respectively. The assumptions show their validity, since brain activity is consistent at both levels compared to previous studies and the time-frequency description correspond to the real frequency values of the EEG bands. The oscillatory responses are stable for each frequency at both levels and as the frequency decreases the stability of the model increases, showing a limit cycle. onclusions
A re-deduction of previous results was made, which had some errors, unifying the notation and many formulas were derived again looking for an algorithmic form that would optimize the calculation. Oscillatory activity occurs when there are interactions between excitatory and inhibitory neuron masses. With the new consideration of delays the temporal dimension was extended, representing the interaction between neurons through a connectivity tensor via specific connection weight and time delays. This new time-delay NMM can simulate several different types of EEG activities since kinetics information was considered, which influence on the oscillatory frequency. An expansion of the previous model was done, considering new structural levels (unit, population and system levels) and different types of human brain rhythms (from gamma to delta). This research makes available a toolbox to move toward more realistic large-scale brain modeling. It allows a high-performance simulation of the neural dynamic behavior according to connectivity information. It enables fast and efficient implementation of those basic core of calculations, due the numerical integration was carried out via algebraic equations using LLM and integrates symbolic and numerical analysis to provide a new way of introducing delays into NMM.
References [1] J. M. Bower,D. Beeman,and M. Hucka. The GENESIS Simulation System California Institute of Technology Pasadena , CA 91125 University of Colorado Boulder , CO 80309 California Institute of Technology [J].2003, (January).
Figure 6: System level simulations. Stable behavior was found also mixing cortical units with different characteristic frequencies.
2] P. A. Valdes-sosa,J. M. Sanchez-bornot,R. C. Sotero,et al. Model Driven EEG / fMRI Fusion of Brain Oscillations [J].2009, 2721:2701 – – –
21. [4] A. Spiegler,S. J. Kiebel,F. M. Atay,et al. Bifurcation analysis of neural mass models: Impact of extrinsic inputs and dendritic time constants [J].NeuroImage, 2010, 52(3):1041 – – – –
18. [9] B. H. Jansen,G. Zouridakis,and M. E. Brandt. A neurophysiologically-based mathematical model of flash visual evoked potentials [J].Biological Cybernetics, 1993, 68(3):275 – F. Grimbert and O. Faugeras. Analysis of Jansen’ s model of a single cortical column [J].Current Biology, 2006, 14(1):34. [11] H. Kwakernaak and M. Sebek. Polynomial J-spectral factorization [J].IEEE Transactions on Automatic Control, 1994, 39(2):315 – – – – – – –
14. [19] I. Bojak and D. T. J. Liley. Axonal velocity distributions in neural field equations [J].PLoS Computational Biology, 2010, 6(1). [20] K. J. Friston,L. Harrison,and W. Penny. Dynamic causal modelling [J].2003, 19:1273 – – –
28. [24] R. Moran,D. A. Pinotsis,and K. Friston. Neural masses and fields in dynamic causal modelling [J].Frontiers in Computational Neuroscience, 2013, 7(APR 2013):1 –
12. [25] A. Spiegler and V. Jirsa. Systematic approximations of neural fields through networks of neural masses in the virtual brain [J].NeuroImage, 2013, 83:704 – ppendix This section collects several equations and mathematical formulations of the LLM due to their length. ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , , ,0 0 10 0 0 l l l n l nn n nl l n l l n l l n p t p tp t p t f p th + − = x p x x x f ff J JJ (1) ( ) ( ) ( ) ( ) ( )
00 12 20 0 0 10 0 0 0 l l l l l l ll l l l l l l l xH p H z dS z t H p S z t y x +− − + − − = f J (2) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) l ll ll l l ll ll h hll h hl lh h h hl l l l ll ll h hl lH hl l ll ll t l l l e h e he e hH e h e H e h e hhH e h eH e h hp t S z tp t z t dS z t − −− −− − − −−− + − = − − − − + − + − + = − + − + = + ABE