Emergence of multifrequency activity in a laminar neural mass model

Emergence of multifrequency activity in a laminar neural mass model | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Emergence of multifrequency activity in a laminar neural mass model View ORCID Profile Raul de Palma Aristides , View ORCID Profile Pau Clusella , View ORCID Profile Roser Sanchez-Todo , View ORCID Profile Giulio Ruffini , View ORCID Profile Jordi Garcia-Ojalvo doi: https://doi.org/10.1101/2025.06.16.659878 Raul de Palma Aristides 1 Department of Medicine and Life Sciences, Universitat Pompeu Fabra , Barcelona, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Raul de Palma Aristides For correspondence: rparistides{at}gmail.com Pau Clusella 2 Department of Mathematics, Universitat Politécnica de Catalunya , Manresa, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Pau Clusella Roser Sanchez-Todo 3 Center of Brain and Cognition, Universitat Pompeu Fabra , Barcelona, Spain 4 Brain Modeling Department , Neuroelectrics, Barcelona, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Roser Sanchez-Todo Giulio Ruffini 4 Brain Modeling Department , Neuroelectrics, Barcelona, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Giulio Ruffini Jordi Garcia-Ojalvo 1 Department of Medicine and Life Sciences, Universitat Pompeu Fabra , Barcelona, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Jordi Garcia-Ojalvo Abstract Full Text Info/History Metrics Data/Code Preview PDF Abstract Neural mass models (NMMs) aim to capture the principles underlying mesoscopic neural activity representing the average behavior of large neural populations in the brain. Recently, a biophysically grounded laminar NMM (LaNMM) has been proposed, capable of generating coupled slow and fast oscillations resulting from interactions between different cortical layers. This concurrent oscillatory activity provides a mechanistic framework for studying information processing mechanisms and various disease-related oscillatory dysfunctions. We show that this model can exhibit periodic, quasiperiodic, and chaotic oscillations. Additionally we demonstrate, through bifurcation analysis and numerical simulations, the emergence of rhythmic activity and various frequency couplings in the model, including delta-gamma, theta-gamma, and alpha-gamma couplings. We also examine how alterations linked with Alzheimer’s disease impair the model’s ability to display multifrequency activity. Furthermore, we show that the model remains robust when coupled to another neural mass. Together, our results offer a dynamical systems perspective of the laminar NMM model, thereby providing a foundation for future modeling studies and investigations into cognitive processes that depend on cross-frequency coupling. Author Summary Understanding how the brain generates and coordinates rhythms across different layers is essential for uncovering the mechanisms underlying perception, memory, and cognition. In this work, we analyze a previously developed model of mesoscopic brain activity that simulates the layered structure of the cortex and its ability to produce coupled slow and fast neural oscillations. Using tools from dynamical systems theory, we reveal how the model gives rise to a rich repertoire of dynamical patterns—including periodic, quasiperiodic, and chaotic activity—through the coexistence of multiple oscillatory modes. We also investigate how pathological changes, such as those linked to Alzheimer’s disease, alter the model’s dynamics and impair its capacity to sustain complex cross-frequency interactions. Finally, we show that the model remains stable when connected to another brain region, highlighting its robustness. Our findings provide a deeper understanding of how multifrequency neural rhythms may emerge, how they might break down in disease, and how this modeling framework can inform both future theoretical studies and the development of new brain models. 1 Introduction Oscillations are ubiquitous in the brain [ 1 ]. These rhythmic or repetitive patterns of electrical activity emerge from the firing of individual neurons, reflecting the synchronized activity of neuronal ensembles, and are therefore the result of the interplay between different scales: microscopic (single neurons), mesoscopic (local groups of neurons), and macroscopic (across brain regions) [ 2 , 3 ]. The frequencies in which mammalian neuronal networks oscillate span four orders of magnitude, ranging from 0.05 Hz to 500 Hz, with this diversity arising from neuron and synapse properties, their interactions, and the circuit motifs that they form [ 4 ]. Neural mass models (NMMs) can capture the mesoscopic dynamics of ensembles of neurons as nonlinear oscillators [ 5 – 9 ]. In this formalism, a population is typically characterized by its average firing rate, membrane potential, as well as the characteristics of synapses through which it connects to other populations (time constant and amplitude of post-synaptic potentials, or PSPs). Although limited, these models have been instrumental in advancing our understanding of brain rhythms, successfully replicating brain activity patterns observed through electrophysiological and neuroimaging techniques [ 10 – 23 ]. In this paper we focus on the emergence of multifrequency coupling in the recently introduced laminar NMM (LaNMM) [ 24 ], an extended NMM that includes the layered structure of the cortex and describes the frequency specificity and coupling between layers. Such a phenomenon appears to be crucial for brain activity, and plays a role in cognitive functions [ 25 – 32 ]. Recent work with the LaNMM has shown how this model can be used to implement gating and predictive coding elements [ 33 ] and cooperation/competition across brain regions [ 34 ]. Frequency specificity is observed across cortical layers, with different frequency bands contributing to distinct aspects of neural processing and communication [ 26 , 35 , 36 ]. For example, during a memory task, macaques exhibit high-frequency activity in superficial layers and low-frequency activity in deep layers. Moreover, the phase and amplitude of deep-layer low-frequency oscillations influence the phase and amplitude of superficial gamma dynamics, highlighting the intricate interplay between frequencies across cortical layers [ 36 ]. From the dynamical point of view, this is the result of the interaction of different nonlinear oscillators, each one with a different natural frequency. In turn, this frequency depends on the intrinsic properties of each neuronal population, such as the synaptic characteristics of neurons that make up the population and the connectivity between them. Hence, when we consider the interaction between different populations, we have to take into account the tendency of each population to oscillate in its preferred frequency. From this interplay, many interesting and complex phenomena emerge, including, multifrequency coexistence, chaos and cross-frequency coupling [ 37 – 42 ]. The LaNMM combines two well-known NMMs: the Jansen and Rit’s (JR) model [ 7 , 8 , 43 ] and a model implementing the pyramidal-interneuron network gamma (PING) mechanism [ 5 , 44 – 47 ]. The first model, consisting of a population of pyramidal neurons coupled to two other populations, is responsible for sustaining lower frequencies and is associated with deeper cortical layers. The latter model includes a pyramidal and an interneuron population, and generates faster rhythms in the gamma band. Finally, inputs to the two pyramidal populations are added to simulate signals from other cortical regions. By integrating these elements, the LaNMM elegantly displays cortical oscillations across different frequency bands [ 22 , 24 ], and has been used to model phenomena including effects of serotonergic psychedelics [ 48 ], oscillatory changes in Alzheimer’s disease [ 33 , 49 ], predictive coding mechanisms [ 33 ], and cooperative vs. competitive interactions [ 34 ]. In particular, the last two studies build on the cross-frequency coupling characteristics of the LaNMM. Here, we analyze the LaNMM’s ability to display multifrequency dynamics as a function of the (constant) inputs of its pyramidal populations. Specifically, we conduct a bifurcation analysis of the model with respect to changes in the input parameters. We show that in addition to the alpha-gamma coupling described in [ 22 , 24 , 33 , 34 ], the model is also capable of generating delta-gamma and theta-gamma coupling, which, to our knowledge, has not been previously reported for other NMMs. Our analysis reveals that the parameter space of external inputs contains a large region of quasiperiodicity, where multifrequency coupling emerges due to the coexistence of two unstable limit cycles associated with the two components of the model: the JR and the PING. Additionally, we examine how alterations associated with Alzheimer’s disease impair this capability, by analyzing their impact on the bifurcation diagram and how they reshape the regions where the system exhibits multifrequency activity. Furthermore, we demonstrate that the model remains robust when coupled to other neural masses. 2 Results 2.1 Laminar Neural Mass Model Neural mass models are lumped mathematical descriptions of neuronal activity that aim to mimic the behavior of populations of neurons with a small number of differential equations. Our main object of study is the LaNMM, first proposed in [ 22 , 24 ], which consists of five interconnected populations of neurons, each representing a different cell type: Pyramidal cells located at deep cortex layers ( P 1 ). Pyramidal cells located at superficial layers ( P 2 ). Fast inhibitory parvalbumin-positive cells ( PV ). Slow inhibitory somatostatin-expressing cells ( SST ). Excitatory spiny stellate cells ( SS ). Additionally, we consider two generic external inputs targeting P 1 and P 2 that represent incoming inputs from other neuronal populations. Figure 1 illustrates the model and the connectivity between these populations. Excitatory synapses, mediated by AMPA neurotransmitters, are represented as triangles, while inhibitory synapses, mediated by GABA neurotransmitters, are shown as circles. Download figure Open in new tab Figure 1. Illustration of the neuronal populations of the LaNMM and the connectivity between them. The top two blocks correspond to the PING model, and the bottom three blocks to Jansen and Rit’s model. Rounded shapes represent interneuron populations, while the rest represent excitatory populations. External inputs targeting P 1 and P 2 are represented by φ e 1 and φ e 2 , respectively. The dynamics of each population and their interactions are modeled following heuristic principles proposed by seminal studies in the field [ 5 , 7 , 50 , 51 ]. First, the average axonal pulses originated in one population (with average firing rate r ( t )) produce an average postsynaptic potential (PSP) y ( t ) via a linear convolution [ 51 ]: The kernel h reads: were parameters A, a vary depending on the effect of the modeled neurotransmitter on the postsynaptic population. For instance, while excitatory neurotransmitters cause excitatory postsynaptic potentials (EPSPs), which depolarize the postsynaptic population, inhibitory neurotransmitters cause inhibitory postsynaptic potentials (IPSPs), which have the opposite effect. Specifically, the amplitude and timescale of the postsynaptic potentials depend on the parameters A (mV), which represents the synaptic gain, and a (s −1 ), which is the time constant of the average PSP. Figure 2(a) shows the three different PSPs captured by Eq. (2) for individual populations (i.e. in the absence of coupling), one for excitatory synapses (mediated by AMPA neurotransmitter, A , blue line) and the other two for inhibitory synapses (mediated by slow and fast GABA neurotransmitter, G s and G f , green and yellow lines respectively). Download figure Open in new tab Figure 2. The three types of h kernels present in the LaNMM lead to different post-synaptic potentials. The PSP amplitudes and decay time are governed by the parameters A and a , respectively. The right hand side displays the frequency response of the filters. In Fig. 2 , the strongest PSPs are generated by GABAergic interneurons, with PV neurons (with GABA f ) responding the fastest and SST neurons (with GABA s ) the slowest. The IP-SPs have a larger amplitude than the EPSPs, due to inhibitory neurons establishing synapses closer to the cell body of postsynaptic cells [ 7 , 52 ]. The frequency response associated with each h-block is shown in Fig. 2 (b), and is modeled using the expression , which describes how each synapse amplifies or attenuates signals across different frequencies. AMPA is the least responsive of all synapses, with a weaker response at higher frequencies. In contrast, fast GABA shows more responsiveness at higher frequencies, i.e. it is more sensitive to fast oscillations and rapid signals, which is typical for fast inhibitory synapses. A second key ingredient of the modeling is a nonlinear function σ ( v ) that relates the average membrane potential ( v ) to an average firing rate [ 53 ]: where φ 0 is the half maximal firing rate of the targeted population, v 0 is the value of the potential when φ 0 is achieved, and r determines the steepness of the sigmoid at the threshold ( v 0 , φ 0 ). Combining these concepts, the average membrane potential y ( t ) of a neuronal population (modulated by interactions with another population through the average membrane potential of the former v ( t )) is governed by the following second-order differential equation [ 54 ] that implements the linear convolution described in Eq. 1 : Using Eq.(4) and considering the connectivity scheme described in Fig. 1 , the system of equations governing the model is given by where the description and values of the parameters are given in Table 1 . In this paper we will pay particular attention to the effects that the external inputs φ e 1 and φ e 2 (acting upon the pyramidal populations P 1 and P 2 , respectively) have on the dynamics of the system. View this table: View inline View popup Download powerpoint Table 1: Description of model parameters and their standard values, as used in [ 22 ]. In the next section we present a detailed analysis of the dynamics and bifurcations of system (5). 2.2 Baseline dynamics of the LaNMM In this section, we explore the dynamical properties of the system described by Eq. (5). We begin by analyzing the model’s response to external inputs. Specifically, we simulate the system starting from a fixed point (with φ e 1 = φ e 2 = 0) and apply a pulse simultaneously to both P 1 and P 2 at t = 0. Figure 3 illustrates the time evolution of the membrane potentials of the model’s two main neural populations. These observables are given by [ 22 ]: Download figure Open in new tab Figure 3. Response of the model for a pulse of 1 ms delivered to P 1 and P 2 simultaneously. Parameters values as in Table 1 , with φ e 1 = 0 and φ e 2 = 0. For better visualization, we center the signals at v = 0 by removing their DC components, to enhance the comparison between responses. We observe differences not only in amplitude but also in the timescale of the response. P 2 exhibits a stronger response, characterized by a fast initial peak followed by damped oscillations back to the steady state. In contrast, P 1 has a significantly slower response with a smaller amplitude. These differences in response dynamics align with the intrinsic properties of these populations: P 2 has higher excitability (lower v 0 ) and faster intrinsic time scales due to its coupling with PV cells. The slower recovery of P 1 is influenced by SST cells, which exhibit the slowest response, as shown in Fig. 2 . On the other hand, the network responses shown in Fig. 3 differ significantly from those observed in Fig. 2 , which is to be expected since populations with distinct amplitudes and time constants now interact. We will revisit this response when analyzing the stability of the fixed point of the system in the next section. As previously mentioned, a key feature of this neural mass model is its ability to sustain distinct oscillatory frequencies across different populations, particularly in the alpha and gamma bands. This behavior is illustrated in Fig. 4 , where the time evolution of the membrane potentials of populations P 1 and P 2 are shown in panels (a) and (b), respectively, for φ e 1 = 200 and φ e 2 = 90 as in the model in [ 22 ]. While both populations exhibit rhythmic activity, their temporal dynamics differ significantly, with P 1 oscillating at a lower frequency, around 10 Hz (alpha range), compared to the faster oscillations of P 2 , around 40 Hz (gamma range). Download figure Open in new tab Figure 4. Time evolution of the membrane potentials of P 1 (a) and P 2 (b) with the respective power spectra shown in panel (c). Parameters values as in Table 1 with φ e 1 = 200 φ e 2 = 90. To further characterize these differences, panel (c) of Fig. 4 displays the power spectra of each signal. The two curves are plotted using the same colors as the corresponding time series in panels (a) and (b). This analysis confirms that while both populations share a strong α component around 10 Hz, P 2 exhibits pronounced γ -band activity at around 40 Hz. The shared α rhythm enables cross-frequency coupling, where γ oscillations in P 2 may be modulated by the phase of the α rhythm. These features were explicitly introduced into the model to ensure a representation of experimental results [ 22 ]. 2.3 Fixed points and bifurcations We now look for the fixed points of the system (5). To solve this system and find the fixed points, we used the continuation and analysis software AUTO-07p [ 55 ]. Figures 5(a) and (b) show a one-parameter bifurcation diagram obtained by varying the external input of the pyramidal population P 1 , φ e 1 , while keeping φ e 2 = 0, in terms of the observables (a) and (b). Download figure Open in new tab Figure 5. Bifurcation diagram of the system showing the fixed points and the amplitude of oscillatory solutions for v P 1 and v P 2 in panels (a) and (b), with φ e 2 = 0. Stable (unstable) fixed points are shown in dark (light) grey. Increasing the external input ( φ e 1 leads to the emergence of oscillatory activity through a SNIC bifurcation (see text). Further increases in φ e 1 result in a Fold of Limit Cycle (FLC) bifurcation, followed by a sequence of supercritical Hopf bifurcations (HB + ). Different oscillatory regimes are marked by asterisks ( * ) and are plotted in panels (c-d), (e-f) and (g-h). Although the three oscillatory regimes are periodic, each one exhibits a different main frequency: in (c-d) θ ≈ 4 Hz, in (e-f) α ≈ 10 Hz, and in (g-h) γ ≈ 40 Hz. Other parameter values are given in Table 1 . Stable (unstable) fixed points are shown in dark (light) grey, together with the amplitude of stable (unstable) periodic solutions in dark (light) colors. The local stability of periodic solutions is given by the Floquet multipliers, which determine if perturbations to such solutions grow or decay. The bifurcations related to the emergence of oscillatory dynamics are marked by vertical lines in panels (a) and (b). These figures show that the system has one stable fixed point for small inputs to P 1 . This stable fixed point, located on the lower branch, has three distinct pairs of complex eigenvalues and two real eigenvalues. This implies that when perturbed, the system’s response will be a combination of oscillatory behavior and exponential decay toward the steady state, in line with the results in Fig. 3 . Increasing φ e 1 leads the system through a saddle-node on an invariant cycle (SNIC) bifurcation at φ e 1 ≈ 107, as two fixed points collapse and a stable invariant cycle appears. Notice that this bifurcation is also encountered in the single Jansen’s model [ 43 ]. The amplitude of the oscillations decreases and the frequency increases after the system passes through two fold (saddle-node) bifurcations of limit-cycles (FLC 1,2 ) at φ e 1 ≈ 159, which are too close to each other to be seen separately in the diagram. This oscillatory activity vanishes at a supercritical Hopf bifurcation at φ e 1 ≈ 367. Increasing φ e 1 leads to a second supercritical Hopf bifurcation at φ e 1 ≈ 457 and to a second invariant cycle, which has a higher frequency a much smaller amplitude (in P 1 ) than the first one. In what follows, we analyze the oscillation dynamics of the model part by part. We divide Figs. 5 (a) and (b) into three different regions marked with * ’s, each one associated with different oscillatory behaviors exhibited by the system. Each one of them is shown in panels (c-g) for v P 1 and (d-h) for v P 2 . The activity in all three regions is periodic and stable, as confirmed by the corresponding Floquet multipliers obtained with AUTO-07p [ 55 ]. In the region labeled as ( * ) in Fig. 5(a) , between the SNIC and the FLC bifurcations, the system displays theta rhythmic activity with frequency ≈ 4 Hz. It is important to notice that although the main frequency of both signals is the same, v P 2 exhibits a fast oscillation on top of the main one. In the region between FLC and the , labeled as ( ** ) in Fig. 5(a) , the system exhibits alpha rhythmic activity with frequency ≈ 10 Hz, and similarly to the region before, v P 2 shows a fast oscillation for v P 2 superimposed on the slower one. The results described above reveal that not only exhibits a fast oscillation, but this oscillation is locked with the phase of . This can be understood by considering the two parts of the model, the JR and the PING, as two coupled oscillators. In regions * and ** , the JR drives the PING activity. This driving is modulated by the function ( Eq. (3) ), which increases considerably as approaches v 0 = 1. Lastly, gamma rhythmic activity is observed in the region labeled as * * * , after the bifurcation, with frequency ≈ 40 Hz. This can be explained by analyzing the dynamics of the single Jansen’s model [ 43 ]. For this input value, the system does not oscillate; instead, it settles into a stable fixed point. In contrast, the PING model exhibits oscillations. As a result, the JR model is being driven by the PING model, but its inherent dynamics suppress the oscillations, trying to pull the system toward its fixed point and hence leading to small amplitudes in the P 1 oscillations. Summarizing, for φ e 2 = 0 the system does not exhibit simultaneously multifrequency activity, in the sense that both v P 1 and v P 2 share a single dominant frequency. That is, only one main frequency is observed across both populations. In the following, we analyze how increasing φ P 2 alters this scenario and can lead to the emergence of distinct dominant frequencies for each population. The effect of a non-zero input in the P 2 population is shown in Figure 6 , which presents different bifurcation diagrams with increasing φ e 2 from left to right, using the same color scheme as Fig. 5 . As we see in the first column (panels a,d), increasing the value of φ e 2 to 20 makes little difference in the position of the first bifurcations (SNIC, FLC, and ). However, the bifurcation point has moved significantly to the left. As we increase φ e 2 = 40, we notice a change in the sequence of bifurcations. Now, after the FLC bifurcations, the system goes through a subcritical Hopf bifurcation (the former bifurcation becomes subcritical when it crosses ), which leads to an unstable periodic solution. Download figure Open in new tab Figure 6. Effect of increasing φ e 2 on the bifurcation diagram. Bifurcation diagrams obtained by varying φ e 1 for different values of φ e 2 . Increasing φ e 2 (left to right) lowers the threshold for the second Hopf bifurcation, , which appears at smaller values of φe 1. Additionally, it also induces a torus bifurcation (TR). Other parameter values are given in Table 1 . Notice that from this point on, until the system reaches , there is an overlap between former regions ** and *** , sustaining alpha and gamma rhythms, respectively. Between these two bifurcations, we notice another one, namely a torus (or secondary Hopf) bifurcation, TR (yellow dashed line in panel (c-f)). Before the TR bifurcation, the unstable inner limit cycle repels the trajectories that are attracted by the outer stable limit cycle (panel (c)). The crossing of the TR bifurcation leads to the loss of the stability of the periodic solution, leading to the emergence of the quasiperiodic behavior. In this regime, the trajectories are attracted to a two-dimensional torus and repelled from the inner limit cycle. For φ e 2 = 90, the scenario remains qualitatively similar, with both HB − 2 and TR shifting significantly to the left, which enables a larger region of overlap between the different limit cycles. Notice that the values used in Fig. 4 lie on the region slightly after the TR bifurcation. The two-parameter bifurcation diagram ( φ e 1 × φ e 2 ) shown in Fig. 7 gives us a better view of how the bifurcation points change for different parameter values. We see that except for the TR bifurcation, the other bifurcations are not significantly affected by the increase of φ e 2 , in agreement with the results of Fig. 6 . To better understand how the bifurcation structure affects the frequency of the signals and , we numerically solved Eq. (5) varying φ e 1 and φ e 2 , with results in presented in Fig. 7 panels (b-c). By comparing the main frequencies of and (panels (b) and (c), respectively) we verify the presence of regions where both populations display the same frequency, such as before the SNIC and after the bifurcations, and also regions with multifrequency activity. Download figure Open in new tab Figure 7. Two-parameter bifurcation diagram of Eq. (5). The color scheme used for bifurcations is the same as in Fig. 2 . Panels (b) and (c) show regions colored according to the dominant frequencies of (b) and (c) , highlighting the model’s ability to sustain multifrequency activity depending on the combination of φ e 1 and φ e 2 . The area of multifrequency activity is bounded by the SNIC (blue) and (purple) bifurcations. In panel (b), we observe that the dominant frequency of remains robust as φ e 2 increases. Between the SNIC and FLC bifurcations, exhibits oscillations in the delta (0.5–4 Hz) and theta (4–7 Hz) ranges. As we move from the FLC to the bifurcations, transitions to oscillations in the alpha range (8–13 Hz). For (panel c), the scenario differs. As φ e 2 increases, the system undergoes a TR bifurcation, shifting the regime and causing the main frequency to move from the alpha range to the gamma range ( > 30 Hz) between the FLC and bifurcations. Additionally, in the region between the SNIC and FLC bifurcations, increasing φ e 2 induces a change in the rhythmic activity of , transitioning from delta/theta frequencies to gamma frequencies for φ e 2 > 250. Other parameter values are given in Table 1 . Specifically, we identify two regions of multifrequency activity. The first occurs for φ e 2 > 250, bounded by the SNIC and FLC bifurcations. In this region, the system exhibits frequencies ranging from approximately 2.5 to 5.0 Hz in , spanning both δ (0.5–4 Hz) and θ (4–8 Hz) bands, depending on the value of φ e 1 . Concurrently, displays γ band activity. The second region of multifrequency activity arises as φ e 1 increases beyond the FLC bifurcation. Here, transitions into the alpha band until it reaches the bifurcation, while maintains γ frequencies (above 30 Hz) after crossing the TR bifurcation. This results in a second distinct region where the two populations oscillate at separate dominant frequencies. These different dynamics are highlighted in Fig. 8 , where the time evolution of the system corresponding to the ( φ e 1 , φ e 2 ) values indicated by the markers ( × , +, ∘) in Fig. 7 . Panels (a) and (b) of Fig. 8 show that δ activity is displayed by P 1 , while γ activity is observed in P 2 . Increasing φ e 1 results in an increase of the frequency of to the θ range, and while this component is also observed in the signal of , the main frequency for this populations remains on the γ range, as shown in panels (c) and (d). Lastly, by increasing φ e 1 further, we show the coupling between α and γ frequencies in panels (e) and (f). Interestingly, we notice that the amplitude of which displays γ rhythm is modulated by the cycle of in all three cases ( δ, θ , and α rhythms). Download figure Open in new tab Figure 8. The model exhibits different frequency couplings depending on φ e 1 and φ e 2 . The time evolution of and is shown for the parameter values highlighted in Fig. (7) , illustrating: (a-b) coupling between δ and γ frequencies, (c-d) coupling between θ and γ frequencies, and (e-f) coupling between α and γ frequencies. All parameter values are specified in Table 1 . We further validate our results by calculating the two largest Lyapunov exponents (LEs) (see Methods), λ 1 and λ 2 , as a function of the external inputs φ e 1 and φ e 2 . The results are shown in Fig. 9 , color-coded according to the values of ( λ 1 , λ 2 ), overlapped with the two-parameter bifurcation diagram. Together, these two values provide insight into the system’s dynamics: The presence of a positive LEs signifies chaos, while negative LEs indicate steadystate dynamics. If one LE is zero, the system exhibits periodic behavior. In contrast, if the system has two zero LEs, it evolves on a two-dimensional invariant torus, indicating quasiperiodicity. As expected, both periodic and quasiperiodic dynamics are observed, with the quasiperiodic regions emerging and being enclosed by torus bifurcations. Interestingly, the system also displays chaotic dynamics in the vicinity of a complex region between FLC bifurcations for high values of φ e 2 . Such a transition from quasiperiodic to chaotic dynamics is also known as torus breakdown . The presence of FLC bifurcations in this region implies the emergence and annihilation of an unstable limit cycle, which can create conditions favorable for chaotic dynamics. Download figure Open in new tab Figure 9. Two-parameter bifurcation diagram of Eq. (5) with regions colored according to the two largest Lyapunov exponents as a function of φ e 1 and φ e 2 . Parameters values as in Table 1 . 2.4 PV interneuron dysfunction In this section we investigate the role of PV cells in the multifrequency observed in the previous sections, more specifically, we focus on reducing the coupling from PV to P 2 , i.e. . Recently, this has been proposed a mechanism to model the accumulation of amyloid-beta oligomers (A β O), which is supposed to damage the synaptic function of PV interneurons [ 49 ]. The accumulation of A β Os in the brain correlates with Alzheimer’s disease (AD) progression, as these soluble oligomers, formed by 2-50 monomers, are believed to be the main culprits behind various neurotoxic effects leading to cognitive decline and begin forming years before clinical signs of the disease appear [ 47 , 56 – 59 ]. Figure 10 reproduces the two-parameter bifurcation diagram from Fig. 7 , focusing on the bifurcation that defines the coexistence of limit cycles for different values of . For comparison, we include the case where , representing the baseline or healthy condition. The region where the two limit cycles coexist is highlighted in grey. Decreasing impacts all the bifurcations. For ( Fig. 10(b) ) and φ e 1 = 0, all bifurcations shift slightly to the left, indicating that oscillatory activity occurs at lower levels of external inputs (hyperexcitability). The most affected bifurcation is the one related to gamma oscillations (HB + , grey curve), whose shift significantly enlarges the grey-shaded area, implying that multifrequency activity is easier achieved. Download figure Open in new tab Figure 10. Bifurcation diagrams showing the effect of reducing on limit cycles and gamma activity. Panel (a) shows the baseline condition, while panels (b)–(d) illustrate the progressive reduction of multifrequency activity and the eventual extinction of gamma oscillations as decreases. The region where the two limit cycles coexist is highlighted in grey. Parameters as in Table 1 . Further decreasing (which would corresponding to increasing the cumulative damage caused by A β Os) reduces gamma activity, as shown in Fig. 10(c) . While the SNIC and bifurcations shift further to the left, the bifurcation shifts to the right, decreasing the grey-shaded area and, consequently, reducing the system’s ability to exhibit multifrequency activity. The results described above suggest that as AD progresses, the regime of co-existence of fast and slow oscillations is reduced, with fast oscillations being disrupted the most. This is important because there is a natural variance of mean inputs to each (LaNMM) brain node as determined by dynamics and tractography. Nodes with lower mean input to P1 and P2 will be affected first. This observation aligns with the evidence on reduced gamma power in AD [ 49 ]. Thus, the model supports the idea that early-stage AD is primarily driven by an imbalance between excitation and inhibition driven by PV interneuron dysfunction [ 49 , 56 , 58 , 60 ]. Additionally, further reductions in result in the complete extinction of gamma oscillations in the LaNMM model, as shown in Fig. 10(d) . 2.5 Long-range connectivity So far, we have considered a model with short-range connections, where all possible external inputs due to long-range connections are encoded by the parameters φ e 1 and φ e 2 . In this section, we briefly investigate the impact of long-range connections on the dynamics of the model, particularly on its ability to sustain multifrequency activity. The coupling scheme between columns is based on the model proposed in [ 61 ] and supported by [ 62 – 64 ]. The final model is illustrated in the center panel of Fig. 11 . Download figure Open in new tab Figure 11. Illustration of the two-column model. (a, b) One-parameter bifurcation diagrams for the oscillators of and as a function of φ . (c, d) Schematic illustrations of the two coupled columns. (e, f) One-parameter bifurcation diagrams for and as a function of φ . The colored lines correspond to the respective pyramidal populations. α = 0.6. Two column model This model represents two cortical columns, each consisting of two layers: a superficial layer which is modeled by the PING part of the LaNMM, and a deep layer modeled by the Jansen’s part of the LaNMM. Long-range connections represent the connectivity between the two cortical columns and are divided into feedback and feed-forward projections. Feedback projections originate from the superficial layer in column 1 and target the superficial layer in column 2. Feed-forward projections originate from the deep layer in column 2 and target both superficial and deep layers in column 1. For simplicity, we assume that populations of the same type have identical parameters A and a across columns. We also consider that all four pyramidal populations receive the same external input, φ . The equations governing the final two-column model are provided in the Methods section. Building on the analysis from previous sections, we now examine how external inputs to the pyramidal populations influence the stability of fixed points and the occurrence of bifurcations focusing on the average membrane potentials of the pyramidal populations. The results are presented in panels (a,b) of Fig. 11 for the first column (c) and in panels (e,f) for the second column (d). The bifurcation diagrams corresponding to the lower layer population ( P 1 ) are shown in the bottom panels, while those of the upper layer population ( P 2 ) are displayed in the top panels. For convenience, the populations are labeled as P 11 and P 12 for the first column (corresponding to P 1 and P 2 of the first column), and P 21 and P 22 for the second column (corresponding to P 1 and P 2 of the second column). Although the introduction of coupling between columns increases the complexity of the system compared to the single-column model, leading to richer dynamical behavior, the bifurcation scenario still bears an important resemblance with the one of Fig. 6 : the coexistence of multiple limit cycles. Before the emergence of any limit cycle in the bifurcation diagrams of Fig. 11 , we notice that the system exhibits one additional stable fixed point compared to the original model, leading to bistability. The emergence of the first limit cycle happens at φ ≈ 35 through an bifurcation at φ ≈ 35, which is unstable until it goes through a TR bifurcation around φ ≈ 332. Further increasing φ leads to another HB − near φ ≈ 86, where an unstable limit cycle emerges. This is followed by an HB + , associated with a stable limit cycle, as shown in the inset of Fig. 11(a) . These two limit cycles collide in a fold of cycles bifurcation (FCIC), giving rise to a stable limit cycle. Once again, the system’s main frequency remains below α , while the P 2 populations exhibit a γ component, which allows δ - γ and θ - γ couplings. As φ continues to increase, the system undergoes a series of period-doubling (PD) and FLC bifurcations around φ ≈ 142. Following these transitions, both columns exhibit α - γ coupling. The multifrequency activity ceases when the limit cycle originating from the FCIC bifurcation disappears in an near the TR bifurcation. Beyond this point, both columns oscillate solely at the γ frequency. The asymmetry in the coupling causes the dynamics of the populations in different columns to differ. For instance, compared to the original model, P 11 receives an extra input from P 21 , and two populations projecting into P 11 ( P 12 and SST 11 ) also receives an input from P 21 . Across the range of φ studied, P 11 exhibits a larger amplitude than P 21 , which reflects the larger number of pyramidal neurons being activated in P 11 . For the top layer pyramidal populations, P 12 and P 22 , although the amplitudes are similar, the dynamics for the range between the SNIC and the FLC bifurcations is rather different. The time evolution of the membrane potentials is shown in Fig. 12 . The top row displays the dynamics of and , while the bottom row shows and . As the external input φ increases from φ = 100 (a–b) to φ = 125 (c–d) and φ = 190 (e–f), the activity of the P 1 populations shifts from the delta to the alpha frequency band. A similar trend is observed in and , although with a more pronounced gamma component. This indicates that while the coupling between columns adds complexity to the model, as observed in the bifurcation diagram, the core feature of the model remains robust and is governed by the same mechanism: the coexistence of limit cycles. Download figure Open in new tab Figure 12. Time evolution of membrane potentials in the two-column model. Top row: Dynamics of and as the external input φ increases, illustrating a transition from delta to alpha frequency bands. Bottom row: Corresponding dynamics of and , showing similar trends but with a prominent gamma-band component. (a–b): φ = 100; (c–d): φ = 125; (e–f): φ = 190. 3 Discussion In this paper we investigate the dynamical properties of a laminar neural mass model, designed to capture different dynamics across cortical layers. The model achieves this by combining two neural mass models, each exhibiting oscillatory dynamics in distinct natural frequency bands, meaning that under external perturbations, each component oscillates at different frequencies. As in the original work, we consider that both pyramidal populations are targeted by external inputs φ e 1 and φ e 2 . To analyze the system’s behavior, we first use φ e 1 as the bifurcation parameter while keeping φ e 2 = 0, allowing us to examine its role in triggering oscillatory dynamics. As shown in Fig. 5 , the oscillatory regime can be divided into three regions depending on the value of φ e 1 , with the dominant frequency transitioning from θ to γ . Additionally, for intermediate values of φ e 1 , the oscillatory behavior resembles that observed in the Jansen model, whereas for higher values, the activity accelerates, corresponding to a PING mechanism for gamma activity. We also investigate the role of φ e 2 in the emergence of rhythmic activity. This input directly targets the PING component of the model, and thus we expect it to facilitate the emergence of gamma oscillations. This is precisely what we observe in Fig. 6 : as we increase the value of φ e 2 , the HB bifurcation that gives rise to gamma activity occurs at smaller values of φ e 1 . This clarifies the mechanism behind the multifrequency activity observed in LaNMM: the coexistence of limit cycles, one sustaining low-frequency oscillations and the other supporting high-frequency oscillations. We verify the robustness of this mechanism through a two-parameter bifurcation analysis as a function of φ e 1 and φ e 2 , shown in Fig. 7 , where we observe a large region in the parameter space where the two limit cycles coexist. Of particular interest are the regions delimited by the torus bifurcations (TR), where the system exhibits quasiperiodic dynamics, and both α and γ frequencies have significant power. Within this region, the model is capable of exhibiting multifrequency activity across different frequency ranges, including δ - γ, θ - γ , and α - γ . We conclude our analysis by evaluating the Lyapunov exponents of the model, shown in Fig. 9 . While quasiperiodicity is observed over a large area of the parameter space, we also identify periodic dynamics in regions where both limit cycles coexist. This behavior can be understood by considering the two layers of the model as entrained oscillators. In such scenarios, the oscillators can resonate, a phenomenon that gives rise to Arnold tongues. This is precisely what we observe in the regions of periodic activity. Moreover, we also observe chaos small regions in the parameter space, in line in recent results. We highlight the significance of the region where the system exhibits two limit cycles by modeling PV interneuron dysfunction, recently proposed as a mechanism to represent A β O, a key biomarker of Alzheimer’s disease progression. As PV cell connectivity deteriorates, the system gradually loses its ability to generate multifrequency activity and, in advanced cases, ceases to oscillate altogether. This suggests that as the disease progresses, the coexistence of fast and slow oscillations diminishes, with fast oscillations being the most affected, in line with the results in [ 49 ]. We also extended the model to account for long-range connectivity. By incorporating a feedforward and feedback coupling scheme between columns, we demonstrate that the coexistence of limit cycles persists under these conditions. This is confirmed through a bifurcation analysis of the extended model, as shown in Fig. 11 . As expected, this model exhibits a more complex bifurcation structure; however, more importantly, rhythmic activity emerges in both the bottom and top layers. As a result, multifrequency activity is observed in both columns at different frequency ranges, as illustrated in Fig. 12 . We emphasize that the mechanism of coexisting limit cycles underlies the approach of incorporating different neural populations to account for multiple frequency ranges in neural mass models, a concept that is not new. The strength of the LaNMM model lies in its ability to exhibit diverse frequency couplings using only five neural populations, controlled solely by adjusting the level of external inputs. While previous models have also been capable of generating multiple frequency interactions, they typically required either increasing the number of neural populations or modifying the time constants of specific populations ( a and φ 0 in our case) to access different frequency regimes [ 41 , 65 , 66 ]. This demonstrates the efficiency of LaNMM in capturing complex neural dynamics with a relatively simple yet powerful parametrization. In order to further refine the LaNMM, future work should generalize our formalism towards next-generation neural mass models [ 9 ]. These modeling framework provides additional equations, rigorously derived from microscopic dynamics, for each neural population mean firing-rate and membrane potential, to which the sigmoid function employed here represents a quasi-stationary approximation [ 67 , 68 ]. To the best of our knowledge, multifrequency activity in these models has been studied only for the case of external forcing [ 69 ], and two interacting inhibitory populations with different time scales [ 70 ]. Our analysis highlights the potential of the LaNMM to model various phenomena and different brain regions in future studies. The ability to exhibit multifrequency activity is a prerequisite for a model to capture cross-frequency coupling phenomena such as phase-amplitude coupling (PAC) or amplitude-amplitude coupling (AAC) and, features commonly observed in neural recordings. For instance, the three frequency couplings observed here are also found in brain activity: δ - γ coupling has been proposed as a biomarker of postictal activity [ 71 ] and is linked to dopamine modulation [ 72 ]; θ - γ coupling is observed in the hippocampus of both rodents and humans [ 73 , 74 ]; and α - γ coupling has been associated with cognitive processes [ 36 , 75 , 76 ]. Recent results show that the cross-frequency coupling exhibited by the LaNMM can provide the key mechanisms for predictive coding, gating [ 33 ] and cooperation/competition across brain regions [ 34 ]. Funding Raul de Palma Aristides, Pau Clusella and Jordi Garcia-Ojalvo are supported by the European Commission under European Union’s Horizon 2020 research and innovation programme Grant Number 101017716 (Neurotwin). Jordi Garcia-Ojalvo was also financially supported by the Spanish Ministry of Science and Innovation, the Spanish State Research Agency and FEDER (Project Reference No. PID2021-127311NB-I00), and by the ICREA Academia program. Giulio Ruffini and Roser Sanchez-Todo are funded by the European Commission under European Union’s Horizon 2020 research and innovation programme Grant Number 101017716 (Neurotwin) and European Research Council (ERC Synergy Galvani) under the European Union’s Horizon 2020 research and innovation program Grant Number 855109. 4 Methods 4.1 Single column model Following the diagram shown in Fig. 1 (a) of [ 22 ], the LaNMM is given in the synapse-driven formalism by With Where s is the type of neurotransmitter and i the population in question. Notice that has a dimension of frequency (Hz). The average membrane potentials are given by: Using the the above equations and taking out the repeated synapses, and since all the synapses coming from one population have the same dynamics and the only difference is the coupling strength relative to the targeting population, we get Notice that u 4 , u 5 , u 12 , and u 13 all represent the perturbation caused by P 1 , which will be weighted by the average number of connections between P 1 and the targets, u 4 . Similarly, u 6 = u 9 = u 11 and u 7 = u 10 . Using these identities we can write: For the sake of coherence with the main text we will introduce the change of variables: which gives us: We now introduce a change of variables u i = C j y i , where C j accounts for the connectivity of the incoming synapse. The first equation, represents the process of receiving an average synapse and producing a membrane perturbation from the point of view of the SST population. The incoming average synapse from P 1 is weighted by C 4 and the outgoing average synapse to P 1 is weighted by C 1 . Taking this into account we have For the external inputs, we have Similarly, for s 8 we have .Putting these results together and using the correspondent parameters for A ’s and a ’s we have: where each one of these second-order equations can be rewritten as a system of two first-order equations using . 4.2 Two-column model Following the diagram shown in Fig. 11 (c,d) the equations of the two column model are given by: To balance the influence of additional inputs on the populations, we introduce the parameter α , which controls the strength of these inputs. For instance, in Eq. (46c), the SST 11 population receives excitatory inputs from both P 11 and P 21 , and α modulates the contribution of these inputs. Additionally, we introduce C 14 , which represents the number of synaptic connections between the P 21 and P 11 populations. All numerical simulations and bifurcation diagrams were performed with α = 0.6 and C 14 = 56.25. Note that α = 0.0 corresponds to the case where the two columns are uncoupled. The impact of this parameter on the system’s dynamics will be explored in future work. 4.3 Numerical simulation and frequency analysis Equations (5) were numerically solved using Python’s scipy.solve ivp function [ 77 ], employing a fourth-order Runge-Kutta algorithm, over a 100-second duration after discarding a 10-second transient, with an integration step of dt = 0.001. The frequency analysis for Fig. 4 was performed using Python’s scipy.fft.fft [ 77 ], with 1-second windows, proper windowing and padding, and averaging. For computational efficiency, in Fig. 7 , the entire signal was used with windowing and padding, and the peak frequency was obtained for P 1 and P 2 . 4.4 Bifurcation diagrams All bifurcation diagrams were calculated with AUTO-07p, a publicly available software for continuation and bifurcation problems [ 55 ]. 4.5 Lyapunov exponents Lyapunov exponents are commonly used to assess how perturbations affect the trajectory of a system over time. A system has as many Lyapunov exponents as its dimensions, but by focusing on the two largest, λ 1 and λ 2 , we can classify the system’s dynamics as follows: Fixed points: both exponents negative ( λ 1 , λ 2 < 0). Periodic dynamics: one zero exponent, the rest negative ( λ 1 = 0, λ 2 0, λ 2 ≤ 0). The Lyapunov exponents for Eq. (5) were calculated using the lyapunovspectrum function from the ChaosTools module in Julia, part of the DynamicalSystems library [ 78 , 79 ]. A threshold of | λ k | < 10 −3 was applied to differentiate between zero and non-zero values. Funder Information Declared European Commission under European Union’s Horizon 2020 , 101017716 European Research Council , 855109 Spanish Ministry of Science and Innovation, the Spanish State Research Agency and FEDER , PID2021- 127311NB-I00 Footnotes https://github.com/raulpalmaa/LaNMM References [1]. ↵ György Buzsáki . Rhythms of the Brain. Online edition, Oxford Academic , 1 May 2009 . Accessed 11 Feb. 2025 . New York : Oxford University Press, 2006 . doi: 10.1093/acprof:oso/9780195301069.001.0001 . OpenUrl CrossRef [2]. ↵ B. Hutcheon and Y. Yarom . “ Resonance, oscillation and the intrinsic frequency preferences of neurons ”. In: Trends in Neurosciences 23 . 5 ( May 2000 ), pp. 216 – 222 . doi: 10.1016/s0166-2236(00)01547-2 . OpenUrl CrossRef [3]. ↵ M. A. Wilson and B. L. McNaughton . “ Reactivation of hippocampal ensemble memories during sleep ”. In: Science 265 . 5172 ( July 1994 ), pp. 676 – 679 . doi: 10.1126/science.8036517 . OpenUrl Abstract / FREE Full Text [4]. ↵ György Buzsáki and Andreas Draguhn . “ Neuronal Oscillations in Cortical Networks ”. In: Science 304 . 5679 ( June 2004 ), pp. 1926 – 1929 . doi: 10.1126/science.1099745 . OpenUrl Abstract / FREE Full Text [5]. ↵ H. R. Wilson and J. D. Cowan . “ Excitatory and inhibitory interactions in localized populations of model neurons ”. In: Biophysical Journal 12 . 1 ( 1972 ), pp. 1 – 24 . OpenUrl CrossRef [6]. W. J. Freeman . Mass Action in the Nervous System . Elsevier , 1975 . doi: 10.1016/c2009-0-03145-6 . OpenUrl CrossRef [7]. ↵ B.H. Jansen , G. Zouridakis , and M.E. Brandt . “ A neurophysiologically-based mathematical model of flash visual evoked potentials ”. In: Biological Cybernetics 68 ( 1993 ), pp. 275 – 283 . doi: 10.1007/BF00224863 . url: https://doi.org/10.1007/BF00224863. OpenUrl CrossRef PubMed [8]. ↵ B. H. Jansen and V. G. Rit . “ Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns ”. In: Biological Cybernetics 73 . 4 ( 1995 ), pp. 357 – 366 . OpenUrl CrossRef [9]. ↵ E. Montbrió , D. Pazó , and A. Roxin . “ Macroscopic Description for Networks of Spiking Neurons ”. In: Physical Review X 5 ( 2015 ), p. 021028 . OpenUrl [10]. ↵ György Buzsáki , Costas A. Anastassiou , and Christof Koch . “ The origin of extracellular fields and currents–EEG, ECoG, LFP and spikes ”. In: Nature Reviews Neuroscience 13 . 6 ( May 2012 ), pp. 407 – 420 . doi: 10.1038/nrn3241 . OpenUrl CrossRef PubMed [11]. F. Wendling et al. “ Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition ”. In: European Journal of Neuroscience 15 ( 2002 ), pp. 1499 – 1508 . doi: 10.1046/j.1460-9568.2002.01985.x . OpenUrl CrossRef PubMed Web of Science [12]. Olivier David and Karl J. Friston . “ A neural mass model for MEG/EEG: coupling and neuronal dynamics ”. In: NeuroImage 20 . 3 ( Nov . 2003 ), pp. 1743 – 1755 . doi: 10.1016/j.neuroimage.2003.07.015 . OpenUrl CrossRef PubMed Web of Science [13]. Gustavo Deco et al. “ The Dynamic Brain: From Spiking Neurons to Neural Masses and Cortical Fields ”. In: PLOS Computational Biology 4 . 8 ( 2008 ), e1000092 . doi: 10.1371/journal.pcbi.1000092 . OpenUrl CrossRef PubMed [14]. S. Coombes . “ Large-scale neural dynamics: Simple and complex ”. In: NeuroImage 52 . 3 ( Sept . 2010 ), pp. 731 – 739 . doi: 10.1016/j.neuroimage.2009.11.045 . OpenUrl CrossRef [15]. Antonio J Pons et al. “ Relating structural and functional anomalous connectivity in the aging brain via neural mass modeling ”. In: Neuroimage 52 . 3 ( 2010 ), pp. 848 – 861 . doi: 10.1016/j.neuroimage.2009.12.105 . OpenUrl CrossRef PubMed Web of Science [16]. Dimitris Pinotsis et al. “ Neural masses and fields: modeling the dynamics of brain activity ”. In: Frontiers in Computational Neuroscience 8 ( Nov . 2014 ), p. 149 . doi: 10.3389/fncom.2014.00149 . OpenUrl CrossRef [17]. Maciej Jedynak et al. “ Temporally correlated fluctuations drive epileptiform dynamics ”. In: NeuroImage 146 ( 2017 ), pp. 188 – 196 . doi: 10.1016/j.neuroimage.2016.11.034 . OpenUrl CrossRef [18]. Maciej Jedynak , Antonio J Pons , and Jordi Garcia-Ojalvo . “ Collective excitability in a mesoscopic neuronal model of epileptic activity ”. In: Physical Review E 97 . 1 ( 2018 ), p. 012204 . doi: 10.1103/PhysRevE.97.012204 . OpenUrl CrossRef [19]. E. Köksal Ersöz et al. “ Neural mass modeling of slow-fast dynamics of seizure initiation and abortion ”. In: PLoS Computational Biology 16 . 11 ( 2020 ), e1008430 . doi: 10.1371/journal.pcbi.1008430 . OpenUrl CrossRef [20]. Nicolás Deschle et al. “ On the Validity of Neural Mass Models ”. In: Frontiers in Computational Neuroscience 14 ( Jan . 2021 ). doi: 10.3389/fncom.2020.581040 . OpenUrl CrossRef PubMed [21]. M. Giugliano , M. Negrello , and D. Linaro J. D. Griffiths , S. P. Bastiaens , and N. Kaboodvand . “ Whole-Brain Modelling: Past, Present, and Future ”. In: Computational Modelling of the Brain . Ed. by M. Giugliano , M. Negrello , and D. Linaro . Vol. 1359 . Advances in Experimental Medicine and Biology. Springer , Cham , 2022 . doi: 10.1007/978-3-030-89439-9\_13 . OpenUrl CrossRef [22]. ↵ Roser Sánchez-Todo et al. “ A physical neural mass model framework for the analysis of oscillatory generators from laminar electrophysiological recordings ”. In: NeuroImage 270 ( 2023 ), p. 119938 . OpenUrl CrossRef [23]. ↵ Michael Forrester et al. “ Whole brain functional connectivity: Insights from next generation neural mass modelling incorporating electrical synapses ”. In: PLOS Computational Biology 20 . 12 ( Dec . 2024 ), e1012647 . doi: 10.1371/journal.pcbi.1012647 . OpenUrl CrossRef [24]. ↵ G. Ruffini et al. “ A Biophysically Realistic Laminar Neural Mass Modeling Framework for Transcranial Current Stimulation ”. In: Clinical Neurophysiology 131 . 4 ( Apr . 2020 ), e78 – e79 . issn: 1388-2457 . doi: 10.1016/j.clinph.2019.12.229 . (Visited on 11/16/2024). OpenUrl CrossRef [25]. ↵ Lawrence M. Ward . “ Synchronous neural oscillations and cognitive processes ”. In: Trends in Cognitive Sciences 7 . 12 ( 2003 ), pp. 553 – 559 . doi: 10.1016/j.tics.2003.10.012 . OpenUrl CrossRef PubMed Web of Science [26]. ↵ A. Maier et al. “ Distinct superficial and deep laminar domains of activity in the visual cortex during rest and stimulation ”. In: Frontiers in Systems Neuroscience 4 ( 2010 ), p. 31 . doi: 10.3389/fnsys.2010.00031 . OpenUrl CrossRef PubMed [27]. E. A. Buffalo et al. “ Laminar differences in gamma and alpha coherence in the ventral stream ”. In: Proceedings of the National Academy of Sciences U.S.A . 108 . 27 ( 2011 ), pp. 11262 – 11267 . doi: 10.1073/pnas.1011284108 . OpenUrl Abstract / FREE Full Text [28]. E. Spaak et al. “ Layer-specific entrainment of gamma-band neural activity by the alpha rhythm in monkey visual cortex ”. In: Current Biology 22 . 24 ( 2012 ), pp. 2313 – 2318 . doi: 10.1016/j.cub.2012.10.020 . OpenUrl CrossRef PubMed [29]. T. Ninomiya et al. “ Microcircuitry of agranular frontal cortex: contrasting laminar connectivity between occipital and frontal areas ”. In: Journal of Neurophysiology 113 . 9 ( 2015 ), pp. 3242 – 3255 . OpenUrl CrossRef [30]. André Moraes Bastos et al. “ Visual Areas Exert Feedforward and Feedback Influences through Distinct Frequency Channels ”. In: Neuron 85 . 2 ( Jan . 2015 ), pp. 390 – 401 . doi: 10.1016/j.neuron.2014.12.018 . OpenUrl CrossRef PubMed Web of Science [31]. J. Bonaiuto et al. “ Lamina-specific cortical dynamics in human visual and sensorimotor cortices ”. In: eLife 7 ( 2018 ), e33977 . OpenUrl CrossRef [32]. ↵ K. Johnston et al. “ Alpha oscillations modulate preparatory activity in marmoset area 8Ad ”. In: Journal of Neuroscience 39 . 10 ( 2019 ), pp. 1855 – 1866 . doi: 10.1523/JNEUROSCI.2703-18.2019 . OpenUrl Abstract / FREE Full Text [33]. ↵ Giulio Ruffini et al. “Cross-Frequency Coupling as a Neural Substrate for Predictive Coding: A Laminar Neural Mass Modeling Approach”. In : ( 2025 ). [34]. ↵ Borja Mercadal et al. “ Bridging local and global dynamics: a laminar mechanistic model for cooperative and competitive interactions in the brain ”. In: bioRxiv ( 2025 ). Manuscript in preparation. [35]. ↵ D. Xing et al. “ Laminar analysis of visually evoked activity in the primary visual cortex ”. In: Proceedings of the National Academy of Sciences of the United States of America 109 ( 2012 ), pp. 13871 – 13876 . doi: 10.1073/pnas.1201478109 . OpenUrl Abstract / FREE Full Text [36]. ↵ A. M. Bastos et al. “ Laminar recordings in frontal cortex suggest distinct layers for maintenance and control of working memory ”. In: Proceedings of the National Academy of Sciences of the United States of America 115 ( 2018 ), pp. 1117 – 1122 . doi: 10.1073/pnas.1710323115 . OpenUrl Abstract / FREE Full Text [37]. ↵ C. C. Chen , S. J. Kiebel , and K. J. Friston . “ Dynamic causal modelling of induced responses ”. In: NeuroImage 41 . 4 ( July 2008 ), pp. 1293 – 1312 . doi: 10.1016/j.neuroimage.2008.03.026 . OpenUrl CrossRef PubMed Web of Science [38]. Viktor Jirsa and Viktor Müller. “ Cross-frequency coupling in real and virtual brain networks ”. In: Frontiers in Computational Neuroscience 7 ( 2013 ), p. 78 . OpenUrl [39]. A. C. E. Onslow , M. W. Jones , and R. Bogacz . “ A Canonical Circuit for Generating Phase-Amplitude Coupling ”. In: PLoS ONE 9 . 8 ( 2014 ), e102591 . doi: 10.1371/journal.pone.0102591 . OpenUrl CrossRef PubMed [40]. Alexandre Hyafil et al. “ Neural Cross-Frequency Coupling: Connecting Architectures, Mechanisms, and Functions ”. In: Trends in Neurosciences 38 . 11 ( Nov . 2015 ), pp. 725 – 740 . doi: 10.1016/j.tins.2015.09.001 . OpenUrl CrossRef PubMed [41]. ↵ Roberto C. Sotero . “ Topology, Cross-Frequency, and Same-Frequency Band Interactions Shape the Generation of Phase-Amplitude Coupling in a Neural Mass Model of a Cortical Column ”. In: PLOS Computational Biology 12 . 11 ( Nov . 2016 ), e1005180 . doi: 10.1371/journal.pcbi.1005180 . OpenUrl CrossRef [42]. ↵ V. Painchaud , N. Doyon , and P. Desrosiers . “ Beyond Wilson–Cowan Dynamics: Oscillations and Chaos Without Inhibition ”. In: Biological Cybernetics 116 . 5-6 ( 2022 ), pp. 527 – 543 . doi: 10.1007/s00422-022-00941-w . url: https://doi.org/10.1007/s00422-022-00941-w. OpenUrl CrossRef PubMed [43]. ↵ François Grimbert and Olivier Faugeras . “ Bifurcation Analysis of Jansen’s Neural Mass Model ”. In: Neural Computation 18 . 12 ( 2006 ), pp. 3052 – 3068 . doi: 10.1162/neco.2006.18.12.3052 . OpenUrl CrossRef PubMed Web of Science [44]. ↵ L. S. Leung . “ Nonlinear feedback model of neuronal populations in hippocampal CA1 region ”. In: Journal of Neurophysiology 47 . 5 ( 1982 ), pp. 845 – 868 . doi: 10.1152/jn.1982.47.5.845 . OpenUrl CrossRef PubMed Web of Science [45]. R. D. Traub , J. G. Jefferys , and M. A. Whittington . “ Simulation of Gamma Rhythms in Networks of Interneurons and Pyramidal Cells ”. In: Journal of Computational Neuroscience 4 ( 1997 ), pp. 141 – 150 . doi: 10.1023/A:1008839312043 . OpenUrl CrossRef PubMed Web of Science [46]. Behnam Molaee-Ardekani et al. “ Computational modeling of high-frequency oscillations at the onset of neocortical partial seizures: From ‘altered structure’ to ‘dysfunction’ ”. In: NeuroImage 52 . 3 ( 2010 ), pp. 1109 – 1122 . doi: 10.1016/j.neuroimage.2009.12.049 . OpenUrl CrossRef PubMed Web of Science [47]. ↵ J. Palop and L. Mucke . “ Network abnormalities and interneuron dysfunction in Alzheimer disease ”. In: Nature Reviews Neuroscience 17 ( 2016 ), pp. 777 – 792 . doi: 10.1038/nrn.2016.141 . OpenUrl CrossRef PubMed [48]. ↵ Jan C. Gendra et al. Restoring Oscillatory Dynamics in Alzheimer’s Disease: A Laminar Whole-Brain Model of Serotonergic Psychedelic Effects . Dec . 2024 . doi: 10.1101/2024.12.15.628565 . (Visited on 12/17/2024). OpenUrl Abstract / FREE Full Text [49]. ↵ Roser Sanchez-Todo et al. “ Fast Interneuron Dysfunction in Laminar Neural Mass Model Reproduces Alzheimer’s Oscillatory Biomarkers ”. In: biorxiv ( 2025 ). [50]. ↵ F.H. Lopes da Silva , A. Hoeks , H. Smits , et al. “ Model of brain rhythmic activity ”. In: Kybernetik 15 ( 1974 ), pp. 27 – 37 . doi: 10.1007/BF00270757 . OpenUrl CrossRef PubMed Web of Science [51]. ↵ A. van Rotterdam et al. “ A model of the spatial-temporal characteristics of the alpha rhythm ”. In: Bulletin of Mathematical Biology 44 ( 1982 ), pp. 283 – 305 . doi: 10.1007/BF02463252 . OpenUrl CrossRef [52]. ↵ Eric R. Kandel et al. Principles of Neural Science . Vol. 4 . New York: McGraw-Hill , 2000 . [53]. ↵ F.H. Lopes da Silva et al. “ Models of neuronal populations: the basic mechanisms of rhythmicity ”. In: Progress in Brain Research 45 ( 1976 ), pp. 281 – 308 . OpenUrl CrossRef [54]. ↵ W.J. Freeman . “ Simulation of chaotic EEG patterns with a dynamic model of the olfactory system ”. In: Biological Cybernetics 56 ( 1987 ), pp. 139 – 150 . doi: 10.1007/BF00317988 . OpenUrl CrossRef PubMed Web of Science [55]. ↵ E. J. Doedel et al. AUTO-07P: Continuation and bifurcation software for ordinary differential equations . Available at: http://indy.cs.concordia.ca/auto/ . 2007 . [56]. ↵ Hyowon Chung et al. “Dissociation of somatostatin and parvalbumin interneurons circuit dysfunctions underlying hippocampal theta and gamma oscillations impaired by amyloid β oligomers in vivo”. en . In: Brain Structure and Function 225 . 3 ( Apr . 2020 ), pp. 935 – 954 . issn: 1863-2661 . doi: 10.1007/s00429-020-02044-3 . url: https://doi.org/10.1007/s00429-020-02044-3 (visited on 12/28/2024). OpenUrl CrossRef PubMed [57]. Kyerl Park et al. “ Optogenetic activation of parvalbumin and somatostatin interneurons selectively restores theta-nested gamma oscillations and oscillation-induced spike timing-dependent long-term potentiation impaired by amyloid β oligomers ”. In: BMC Biology 18 . 1 ( Jan . 2020 ), p. 7 . issn: 1741-7007 . doi: 10.1186/s12915-019-0732-7 . url: https://doi.org/10.1186/s12915-019-0732-7 (visited on 12/28/2024). OpenUrl CrossRef PubMed [58]. ↵ Virginia Garcia-Marin et al. “Diminished perisomatic GABAergic terminals on cortical neurons adjacent to amyloid plaques”. English . In: Frontiers in Neuroanatomy 3 ( Nov . 2009 ). Publisher: Frontiers . issn: 1662-5129 . doi: 10.3389/neuro.05.028.2009 . url: https://www.frontiersin.org/journals/neuroanatomy/articles/10.3389/neuro.05.028.2009/full (visited on 12/31/2024). OpenUrl CrossRef PubMed [59]. ↵ Magdalena Martinez-Losa et al. “Nav1.1-Overexpressing Interneuron Transplants Restore Brain Rhythms and Cognition in a Mouse Model of Alzheimer’s Disease”. eng . In: Neuron 98 . 1 ( Apr . 2018 ), 75 – 89.e5 . issn: 1097-4199 . doi: 10.1016/j.neuron.2018.02.029 . OpenUrl CrossRef PubMed [60]. ↵ Laure Verret et al. “Inhibitory interneuron deficit links altered network activity and cognitive dysfunction in Alzheimer model”. eng . In: Cell 149 . 3 ( Apr . 2012 ), pp. 708 – 721 . issn: 1097-4172 . doi: 10.1016/j.cell.2012.02.046 . OpenUrl CrossRef PubMed Web of Science [61]. ↵ Anna Callau Rabasa . “ Long-range brain connectivity and its role in brain network modeling ”. MA thesis. UPC , 2024 . [62]. ↵ Jorge F. Mejias et al. “ Feedforward and feedback frequency-dependent interactions in a large-scale laminar network of the primate cortex ”. In: Science Advances 2 . 11 ( 2016 ), e1601335 . doi: 10.1126/sciadv.1601335 . OpenUrl FREE Full Text [63]. Jorge Jaramillo , Jorge F. Mejias , and Xiao-Jing Wang . “ Engagement of Pulvino-cortical Feedforward and Feedback Pathways in Cognitive Computations ”. In: Neuron 101 . 2 ( 2019 ), 321 – 336.e9 . doi: 10.1016/j.neuron.2018.11.023 . OpenUrl CrossRef PubMed [64]. ↵ Xiao-Jing Wang . “ Theory of the Multiregional Neocortex: Large-Scale Neural Dynamics and Distributed Cognition ”. In: Annual Review of Neuroscience 45 ( 2022 ), pp. 533 – 560 . doi: 10.1146/annurev-neuro-110920-035434 . OpenUrl CrossRef PubMed [65]. ↵ M. Chehelcheraghi et al. “ A neural mass model of phase—amplitude coupling ”. In: Biological Cybernetics ( 2016 ), pp. 1 – 22 . [66]. ↵ M. Chehelcheraghi et al. “ A neural mass model of cross frequency coupling ”. In: PLoS ONE 12 . 4 ( 2017 ), e0173776 . doi: 10.1371/journal.pone.0173776 . OpenUrl CrossRef [67]. ↵ Federico Devalle , Alex Roxin , and Ernest Montbrió. “ Firing rate equations require a spike synchrony mechanism to correctly describe fast oscillations in inhibitory networks ”. In: PLoS Computational Biology 13 . 12 ( Dec . 2017 ), pp. 1 – 21 . OpenUrl CrossRef [68]. ↵ Pau Clusella et al. “ Comparison between an Exact and a Heuristic Neural Mass Model with Second-Order Synapses ”. In: Biological Cybernetics (Dec. 2022 ). issn: 1432-0770 . doi: 10.1007/s00422-022-00952-7 . OpenUrl CrossRef [69]. ↵ Marco Segneri et al. “ Theta-Nested Gamma Oscillations in Next Generation Neural Mass Models ”. In: Frontiers in Computational Neuroscience 14 ( 2020 ), p. 47 . OpenUrl CrossRef [70]. ↵ Andrea Ceni et al. “ Cross frequency coupling in next generation inhibitory neural mass models ”. In: Chaos: An Interdisciplinary Journal of Nonlinear Science 30 . 5 ( May 2020 ), p. 053121 . issn: 1054-1500 . doi: 10.1063/1.5125216 . eprint: https://pubs.aip.org/aip/cha/article-pdf/doi/10.1063/1.5125216/14627750/053121\_1\_online.pdf . url: https://doi.org/10.1063/1.5125216. OpenUrl CrossRef PubMed [71]. ↵ Vasily Grigorovsky et al. “ Delta-gamma phase-amplitude coupling as a biomarker of postictal generalized EEG suppression ”. In: Brain Communications 2 . 2 ( 2020 ), fcaa182 . doi: 10.1093/braincomms/fcaa182 . OpenUrl CrossRef [72]. ↵ Victoria Andino-Pavlovsky et al. “ Dopamine Modulates Delta-Gamma Phase-Amplitude Coupling in the Prefrontal Cortex of Behaving Rats ”. In: Frontiers in Neural Circuits 11 ( 2017 ), p. 29 . doi: 10.3389/fncir.2017.00029 . OpenUrl CrossRef PubMed [73]. ↵ G. Buzsáki et al. “ Hippocampal network patterns of activity in the mouse ”. In: Neuroscience 116 ( 2003 ), pp. 201 – 211 . doi: 10.1016/S0306-4522(02)00669-3 . OpenUrl CrossRef PubMed Web of Science [74]. ↵ R. T. Canolty et al. “ High gamma power is phase-locked to theta oscillations in human neocortex ”. In: Science 313 . 5793 ( 2006 ), pp. 1626 – 1628 . doi: 10.1126/science.1128115 . OpenUrl Abstract / FREE Full Text [75]. ↵ Bradley Voytek et al. “ Shifts in gamma phase–amplitude coupling frequency from theta to alpha over posterior cortex during visual tasks ”. In: Frontiers in Human Neuroscience 4 ( 2010 ), p. 191 . doi: 10.3389/fnhum.2010.00191 . OpenUrl CrossRef PubMed [76]. ↵ H. Park et al. “ Formation of visual memories controlled by gamma power phase-locked to alpha oscillations ”. In: Scientific Reports 6 ( 2016 ), p. 28092 . doi: 10.1038/srep28092 . OpenUrl CrossRef PubMed [77]. ↵ SciPy 1.0 Contributors . “ SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python ”. In: Nature Methods 17 . 3 ( 2020 ), pp. 261 – 272 . doi: 10.1038/s41592-019-0686-2 . OpenUrl CrossRef PubMed [78]. ↵ George Datseris . “ DynamicalSystems.jl: A Julia software library for chaos and nonlinear dynamics ”. In: Journal of Open Source Software 3 . 23 ( Mar . 2018 ), p. 598 . doi: 10.21105/joss.00598 . url: https://doi.org/10.21105/joss.00598. OpenUrl CrossRef [79]. ↵ George Datseris and Ulrich Parlitz . Nonlinear dynamics: A concise introduction interlaced with code. en . Cham, Switzerland : Springer Nature , 2022 . doi: 10.1007/978-3-030-91032-7 . url: https://doi.org/10.1007/978-3-030-91032-7. OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted June 17, 2025. Download PDF Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Emergence of multifrequency activity in a laminar neural mass model Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Emergence of multifrequency activity in a laminar neural mass model Raul de Palma Aristides , Pau Clusella , Roser Sanchez-Todo , Giulio Ruffini , Jordi Garcia-Ojalvo bioRxiv 2025.06.16.659878; doi: https://doi.org/10.1101/2025.06.16.659878 Share This Article: Copy Citation Tools Emergence of multifrequency activity in a laminar neural mass model Raul de Palma Aristides , Pau Clusella , Roser Sanchez-Todo , Giulio Ruffini , Jordi Garcia-Ojalvo bioRxiv 2025.06.16.659878; doi: https://doi.org/10.1101/2025.06.16.659878 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Neuroscience Subject Areas All Articles Animal Behavior and Cognition (7622) Biochemistry (17648) Bioengineering (13870) Bioinformatics (41880) Biophysics (21423) Cancer Biology (18553) Cell Biology (25458) Clinical Trials (138) Developmental Biology (13364) Ecology (19866) Epidemiology (2067) Evolutionary Biology (24290) Genetics (15589) Genomics (22475) Immunology (17711) Microbiology (40327) Molecular Biology (17145) Neuroscience (88472) Paleontology (666) Pathology (2826) Pharmacology and Toxicology (4815) Physiology (7635) Plant Biology (15114) Scientific Communication and Education (2044) Synthetic Biology (4286) Systems Biology (9815) Zoology (2268)