The WC/JR neural mass models have been widely studied from varying perspectives, as discussed above. Even so, it is worth reiterating that the differential equations in NM models have the form of a driven simple harmonic oscillator (Halliday et al., 1997):

Where y(t) is a time-dependent state variable (e.g., voltage), k is a damping coefficient that measures energy dissipation in the system, τ is the time constant characterizing system dynamics (e.g., decay) and F_ext is an external driving force. In an undamped system (k = 0), the oscillation frequency is ω = 1/τ. Non-linearities in such models enter via F_ext. Equation (1) has been widely used and is well-characterized; in particular, the system's oscillation frequency is Halliday et al. (1997):

When k = 2/τ, Equation (2) takes the form used in NM models, as in Equations (3)–(5). Then, ω = 0 and the system has no characteristic frequency and is a critically damped oscillator. Such a system is a “one shot oscillator” that quickly reverts to its initial state and thus does not sustain oscillations on its own—much like a well-tuned door closer. Its sustained activity only appears when driven by external inputs, such as external noise or stimuli, and by the outputs of linked nodes. Those inputs have a non-linear form, generally sigmoidal in shape, and generate non-linear dynamics. Overall, this is reminiscent of a neuron generating an action potential, followed by a refractory period. The neuron responds to external stimuli but, by itself, does not oscillate. Yet, collections of interacting neurons do. The origins of network oscillations in spiking (Wallace et al., 2011) and other systems have been reviewed (Brunel and Wang, 2003; Wang, 2010).

Casting the NM model in the form of Equation (1) provides intuition to guide the tuning of model parameters. The terms on the left are acceleration, damping force, and restoring force, while those on the right are the driving forces (stimuli and interactions). For stable oscillation, a reasonable balance of these terms is required; otherwise, the oscillation rapidly dies out, as it does for many parameter choices. Realizing the simplicity of Equation (1) invites the further analogy of coupled oscillators and the resulting frequency shifts caused by the coupling strengths. Extensions of Equations (1), (2) interacting oscillators have been well-studied for standard oscillator models, with explicit expressions for the resulting frequency shifts and amplitudes available (Jothimurugan et al., 2016). Such studies provide insight into tuning the frequency of each neural mass in the present study (cf. Supplementary material 2.5).

Each anatomical area in the cluster is a node modeled as three neural assemblies, or neural masses (NM), containing excitatory or inhibitory neurons, following WC (Wilson and Cowan, 1972, 1973) and JR (Jansen and Rit, 1995), originally formulated as models of a cortical column (Mountcastle, 1997; Molnar, 2013). The models used a measure of neuronal activity, such as the fraction of excited neurons or firing rate, and have also been interpreted using voltage, as adopted herein. The two views are equivalent given the monotonic relationship between the firing rate and voltage, noting that the signaling rate plateaus, or saturates, as voltage increases. The dynamics of each NM were originally described by two signal processing blocks that accumulate pulses to generate a post-synaptic potential and transform the voltage to a firing rate (Jansen and Rit, 1995). The resulting differential equation is just that of a critically damped linear oscillator. Its sustained activity only appears when driven by external inputs, from external noise or stimuli, and by the outputs of linked nodes. Those inputs have a non-linear form, generally sigmoidal in shape, and generate non-linear dynamics. The resulting coupled differential equations, when cast in the standard form (derivatives on the left and driving terms on the right), are as follows:

Where y_e and y_i are the potentials associated with the excitatory and inhibitory neural populations, respectively, y′ is its first-time derivative, and y" is the second derivative. Each equation is for a single node k (= 1–6), and the sum in Equation (3) is over linked nodes j. y₀ is the potential of the subpopulation of excitatory neurons that are driven by the difference dy = y_e – y_i. The time constants τ_e and τ_i (~ 10 ms) represent the joint effects of “delays associated with the membrane resistance and distributed delays in the dendritic network” (Jansen and Rit, 1995) in the excitatory and inhibitory neural populations. The units throughout Equations (1)–(3) are V/s², although mV/s² is used in practice. The presence on the left of coefficients 2/τ_e and 1/τ²_e indicates that the oscillator is critically damped, as discussed above. The terms on the right are the non-linear feedback and driving inputs; the latter are applied only to the excitatory population. These internal and external driving forces are what cause the system to produce sustained oscillations, overcoming the critical damping. Casting the model in the form of Equations (1), (3)–(5) aids in the tuning of the parameters. It also makes clear the balance of driving forces, particularly between the excitatory and inhibitory sub-populations.

The input p(t) is the background noise due to local cortical activity. The final term in Equation (1) is the summed inputs to the current node (k) from linked nodes (j), weighted by the adjacency, or connectivity, matrix elements a_kj. Numerous authors have discussed the adjacency matrix and its use in neuroscience; that literature is best captured by reviews by Bullmore and Sporns (2009) and Bassett and Sporns (2017) and textbooks by Sporns (2011). In Equation (3), dy_j = y_e – y_i is the output of each linked node j (in mV) and is the difference in voltage outputs of the excitatory and inhibitory neural subpopulations that are presented to the pyramidal subpopulation in each node. Following Jansen and Rit (1995) and Goodfellow et al. (2012), the model output, dy_j, is taken to be a local field potential for that anatomical area. LFP for the cluster is taken as the average of the six contributions. Signal delays can be introduced in this term by using instead a time-delayed dy_j (t – d_ij/v), where d_ij is the inter-node distance, available from atlas coordinates, and here v is the local signal velocity. Local signals were assumed to travel along axonal or other pathways at 1 m/s (Muller et al., 2014), assuming short-range links may be unmyelinated. The coefficients A and B, originally the maximum post-synaptic potential amplitude, also quantify the feed-forward synaptic strengths, while C₁-C₄ describe the feedback strengths to the excitatory and inhibitory populations (cf. Supplementary material 2.2). The function S (units Hz or s⁻¹) is the voltage to firing rate transformation, described below. Parameter values are discussed below and are taken to be broadly consistent with those used in previous studies (Wilson and Cowan, 1972, 1973; Jansen and Rit, 1995; David and Friston, 2003). Derivation of the average synaptic weights w_av from available data is discussed below. Note that there can be a separate w_av for in and out links, so the function S may vary.

These equations were solved numerically by a 4th order Runga Kutta method (Press et al., 2007), with a timestep of 0.1 ms, using MATLAB codes. Recasting each 2nd order differential equation as two 1st order equations improves numerical stability. Small background noise (Gaussian random) was applied; alternatively, the use of uniform or no noise did not affect the results. Signal delays mean that Equations (3)–(5) become delay differential equations (cf. MATLAB Help Menu). They are of a simple type since delays only appear in the driving terms on the rights of Equations (3)–(5), so they are straightforward to solve. Transients died out within 2–3 s of simulation; stimuli were initiated at 4 s, and the simulation ran for a total of 20 s. Longer duration responses were simulated out to 30 or 40 s as a check. Fast Fourier Transforms of model outputs (voltage dy), using the MATLAB Signal Processing Toolbox, were used to find the spectral density distribution of each NM's output to identify prominent spectral peaks (MATLAB findpeaks function) and to quantify the distribution of power across frequency bands.

The voltage-to-rate transformation function S is usually approximated by the convenient sigmoid function, which can be derived from the distribution over voltage thresholds for firing in a neural population of fixed synaptic weight (Wilson and Cowan, 1972). However, measured firing thresholds vary by only 10s% (Henze and Buszki, 2001; Yu et al., 2008), while connectivity data of link number and link weights can vary by 2–3 orders of magnitude, as with the present data. In that case, a distribution of synaptic weights in a neural population with the same threshold, v₀, is appropriate, as originally noted (Wilson and Cowan, 1972). That derivation (Supplementary material 2.1) yields a new transformation function:

Where q_m is the maximum rate, v₀ is the mid-point (50% firing rate) of the sigmoid curve, w_av is the average synaptic weight in each neural population, and r is the inverse of its standard deviation, and also the steepness of the sigmoid. The function erfc is the complementary error function (Arfken, 1970).¹ This is a steeper curve than the usual sigmoidal form (cf. Supplementary Figure S2).

Many simulations use identical neural masses, in which model parameters are identical for each node. Neural masses are used to model cortical areas or columns, so it is evident from anatomical and other observations that the neural masses need to differ. That can be captured in the internal model parameters described above. An evident variation between anatomical areas is in volume V and number N of neurons in each area, followed by the varying number and weights of the in- and out-links, here captured by w_av. All of those are available from the marmoset data used here. The internal time constants, τ_e and τ_i, are not known but generally fall in the range of 10–25 ms in keeping with known neural time constants (Koch et al., 1996). The values of N, w_{av − in}, and w_{av − out} vary significantly across the six nodes of the cluster studied here, so it is likely that the emergent dynamics of those nodes may vary.

Random noise (zero mean and unit variance) was used to model inputs to an NM from surrounding neural regions, as is usually used. It is also known that wave-like oscillations propagate across the brain, some associated with cognitive tasks (Bhattacharya et al., 2022). Neural field theory provides a detailed description of electromagnetic waves in the brain that manifest as standing waves, or eigenmodes (Robinson et al., 2016), and traveling waves (Gabay et al., 2018). Waves communicating between distant brain regions have been observed in marmoset (Davis et al., 2019). Propagating beta band LFP waves have been shown to facilitate information transfer in motor cortical areas of macaque (Rubino et al., 2007) and to be present in the human motor cortex (Takahashi et al., 2011). Theta and alpha traveling waves have also been observed in the human neocortex, correlated with a memory task (Zhang et al., 2018). In view of this and other evidence, simulated wave-like stimuli were applied to the pFC cluster.

To simulate a wave like stimulus, a signal comprising a constant pulse rate, modulated by a sin wave in one of the four bands, was constructed. The stimulus, being a pulse rate, must be strictly positive (cf. Supplementary Figure S9). The options to achieve this are explained in Supplementary material 2.9. The model was stimulated by pulse train waveforms with amplitude modulated by an oscillatory waveform in each of the theta to gamma bands. The frequency was chosen to match the spectral peaks in the resting state simulation (cf. Table 1). The stimulus was turned on at 4.0 s, well-after transients had settled down, and generally applied either for the full simulation run (20 s). When present, a second stimulus was applied at 10 s. Stimuli were applied to all nodes or to nodes singly, in pairs, etc. When needed, simulation runs were extended to 30 or 40 s if further responses persisted. A pulse function stimulus was also constructed following the original methods (Jansen et al., 1993), using a pulse density function of a monophasic exponential form: t⁷exp(–t/t), with a 4-ms time constant. Its maximum amplitude of 200 Hz matched that of the wave stimulus and delivered 5.2 “unit pulses” (unitless: Hz ^* s) as a proxy measure of energy input. For comparison, a 16-s wave-modulated stimulus typically delivers ~1,000 unit pulses over the duration of the stimulus.

The original NM model, which uses the sigmoid voltage-rate transformation function (Supplementary Equation S2), produces LFP outputs that are combinations of simple sinusoidal functions, as illustrated in Supplementary Figure S5. Some heterogeneity is available via the nodes' internal parameters. Furthermore, mild heterogeneity can be introduced by modifying the sigmoid to include average synaptic weights, w_av (Supplementary Equation S3), yielding similar results. These simple models did not produce dynamic transitions in response to stimuli. The extended NM model, used herein, uses integration over a normal distribution of synaptic weights to yield a steeper, shifted voltage-rate transformation function (Equation 6; Supplementary Equation S4; Supplementary Figure S2) that is more sensitive to w_av and is the focus of the present study. This w_av can be estimated from available linkage data for the marmoset cortex (Supplementary material 2.4; Supplementary Tables S1, S2). Searching parameter combinations produced both simple and complex dynamics. Two (#2, #6 out of 10) trials yielded sustained transitions in the dynamical state following wave-like stimuli of selected nodes in the cluster. The results presented here are for a single set of parameters (trial #6) using the parameters listed in Table 1. A typical output of the model in a resting state, driven only by noise, is shown in Supplementary Figure S7. The LFP for the cluster has spectral power distribution: 63% theta, 14% alpha, 21% beta, and 1.7% gamma bands. A representative Fourier spectrum is shown in Supplementary Figure S8 and is compared to available data, which shows strong beta peaks.

Stimuli were applied to all nodes and selectively to single and pairs of nodes. Node 2 (A32V) showed a distinctive response, presented below, while the other target nodes showed a small (~2 mV) step up in potential with an unchanged waveform. Other nodes, not targeted, showed no response. A gamma (32.1 Hz) stimulus to both A32V and A11 caused the responses shown in Figure 2. The biased full wave and the half wave stimulus caused identical responses, with varying delays in onset. The half wave stimulus (Supplementary Figure S9B) induced a complex, ramped, beating waveform after a 1.6-s transition (Figure 2A), repeating after 6.6 s; while the full wave stimulus induced a simple beat (Figure 2B) after a 2.6 s transition. Note that the full wave stimulus (Supplementary Figure S9A) is effectively oscillating at twice the input frequency (i.e., 64.2 Hz) while still in the gamma band. A check, by applying a double frequency (64.2 Hz) half wave stimulus to both A32V and A11, confirmed a response (not shown) similar to that presented in Figure 2B, but with a longer repeat time of 8.6 s.

A full-wave gamma stimulus applied only to A32V caused a negligible response, while that stimulus only to A11 rapidly (in 0.86 s) induced a different beat pattern (not shown), repeating at 8.6 s. The response of A32V is dominated by the response of the inhibitory subpopulation oscillation, shown in Supplementary Figure S10. All subsequent results are in response to half-wave stimuli (Supplementary Figure S9B). A comparison of responses in A32V to stimuli to both A32V and A11 in the four standard frequency bands is shown in Figure 3. The stimulus frequency in each band is chosen to match the fundamental resonances listed in Table 1.

The A32V waveform transitions to a lower amplitude steady state after varying transition times, with a long delay of 12.5 s for a beta-modulated stimulus (Figure 3C) but only 1.6 s for gamma (Figure 3D). The gamma stimulus induces an additional long period, ramped pattern, repeating after 6.6 s (0.15 Hz), which is discussed below (Section 3.3). The theta stimulus, possibly modeling the multiple inputs of neighboring theta oscillators in the dominant spectral band, induces a long-lasting, low state that eventually (after 17 s) relaxes back to the resting state (cf. Supplementary Figure S11). The alpha stimulus appears to induce a disordered state (Figure 3B) of uncertain interest. Note that the theta stimulus was tuned to the resonance of A32V (5.3 Hz); the other theta resonance, for A46D (5.9 Hz), results in a fast (0.28 s) transition to a multi-beat waveform (not shown), possibly suggesting a different modulating role for A46D. Many tests indicated that consistent results were obtained when stimuli were applied simultaneously to both the in-hub, A32V, and to the in-connector, A11, suggesting a local circuit functioning as an AND logic gate. Network analysis and path tracing (Pailthorpe, 2024) also showed that both these nodes are the ultimate target of many sensory pathways. Stimulus to a single node was tested, e.g., application of the beta wave to A32V only induced a fast (1.8 s) transition to a smaller (by 8 mV) amplitude beat, while beta stimulus to A11 only produced a negligible response (~1 mV change). Generally, the average LFP for the cluster exhibited only small changes. The prominent response was for the local potential at A32V. As a final check, switching off the stimulus wave generally caused the oscillation to revert to its original state.

For the resting state cluster (cf. Supplementary Figures S7, S8), the dominant spectral power is in the theta band (63%) and then in the beta (21%) and alpha (14%) bands. Current pFC reviews (Lundqvist et al., 2018a,b) suggest that transients in beta and gamma bands have a role in working memory readout. Beta traveling waves have been observed to originate from pFC during task performance in macaques (Bhattacharya et al., 2022). To prove this, a beta wave was applied to A32V and A11 at 4 s, and then a second stimulus was applied 6 s later. The reset stimulus was applied to all nodes. Sensitivity was tested by also applying it selectively to single and pairs of nodes. Beta emissions from the cluster were not sought nor observed. The first beta stimulus (Figure 3C) appears to induce a latent switch that, after 12.5 s, transitions to a “low” amplitude state that can be interrupted within 0.1 s by switching off the beta wave at an intermediate time (e.g., 10 s), inducing a return to the original oscillation mode (not shown). That transition also is inhibited by a second, short stimulus, as shown in Figure 4. This transition is affected by a 100-ms gamma burst to the single nodes A32V, A10, or A46D but not by the stimulus of A11 alone.

A 100-ms gamma burst (Figure 4A) applied at 10 s (i.e., before the anticipated transition at 12.6 s) rapidly (~0.1 s) induces a reversion, within 1.0 s, to the initial state (overlaid in gray in Figure 4A). This short stimulus delivers little energy to the oscillators, so it needs to be close to in-phase with the excitatory potential y_e (of A32V) of the responding node to effect the transition; if it is out of phase, a proportionally reduced effect follows. By contrast, a simple pulse function stimulus, given by Equation (3) of Jansen et al. (1993) (reproduced in Section 3.5), interrupts the intermediate oscillation and initiates (within 0.1 s) the transition to the low state (Figure 4B), which then slowly relaxes to the original oscillation. These two stimuli have quite different energy content since a 100-ms gamma burst contains up to 6 positive peaks, each of which is approximately equivalent to a single pulse.

Since the theta band dominates the power spectrum in the resting state, comparison simulations were run: a sustained and a short gamma burst were also applied after the theta-induced transition (Figure 3A). A continuous gamma wave applied at 10 s, after the theta wave induced low state had stabilized, induces a new oscillatory mode that repeated at ~6.7 s, as shown in Figure 5. By contrast, a 100-ms gamma burst had only a small effect on the theta-induced oscillation, even when care was taken to ensure that the stimulus was in phase with the y_e (A32V) oscillation.

Given the interest in gamma bursts, the results were investigated further. The system is a set of physical oscillators, and these transitions are likely to require sufficient energy input to drive them into a new state. Thus, longer gamma bursts and repeated short bursts that were phase-aligned to the positive parts of the theta oscillation were simulated. Recall that the theta-modulated dynamics spontaneously revert to the high amplitude state after 17 s (cf. Supplementary Figure S11). A single 500-ms gamma burst induced that reversion after 14 s, and a double amplitude 500 ms burst reduced that wait time to 12.3 s (figure not shown). Such a longer burst drives the oscillators with, and then against, the original theta-induced driving force. To avoid that, three 100-ms gamma bursts were constructed to phase align with only positive cycles of the theta oscillation that induced the transition after 13.9 s. This illustrates how gamma bursts, properly phase aligned, can deliver more efficient energy transfer and deliver a faster dynamical state transition of the cluster.

A number of structural clusters in the marmoset cortex are described in my earlier study (Pailthorpe, 2024), with the pFC cluster being the most distinctive. It is compact and almost fully linked, with 29 of 30 possible internal links present; one-third of those have strong weights. Other clusters were also investigated amongst the remaining hubs and prominent connector nodes in the cortical network. Two candidates were an auditory group (AuA1, AuCPB, and AuRT) and a visual group (V3, V4, and V6). A motor-somatosensory group [cf. Supplementary Table S2 of Pailthorpe (2024)] was found to be more spatially dispersed and less tightly interlinked; another in-hub, TEO, was linked with the visual group. These clusters were weak, longer range, and less fully interlinked. Their linkage patterns were more tree-like rather than the tight, star-shaped patterns found in the pFC cluster. All these other clusters warrant further investigation.

All six nodes in the cluster received inputs from the visual cortex, mostly from V2 via TPO. Other links, via MST or Opt or direct from V2 or V4, are weaker. The strongest weight pathways from the auditory cortex are from AuA1 and AuCPB via TPO to A32V and A11, consistent with the observations of Dorsal and Ventral auditory pathways (Schreiner and Winer, 2007) and the sensitivity of A32V and A11 to stimuli in the simulations. The somatosensory areas (cf. Section 4.1, above) have direct links to A10 and A11. There are reciprocal links between A10 (the out hub) and the pre-motor area A6DR, which has a strong input to A46D in the pFC cluster. The strong links to A46D, along with its simulated responses, suggest it may have a role in modulating pFC cluster responses.

The present study demonstrates dynamical transitions in localized clusters of cortical areas modeled as simple oscillators. On a larger scale, two studies shed light on bursts, transitions, and synchronization and are accompanied by significant literature. The first was an experimental investigation of spontaneous LFP oscillations in rat somatosensory cortex measured by micro-electrode arrays (Giresch and Plenz, 2008), which observed beta and gamma bursts nested within theta oscillations. The oscillations were coherent between nearby areas, with the synchronization following a power law across distances up to 1.6 mm, and the nested oscillations being organized as avalanches. In the present study, higher frequency, brief stimuli needed to be phase-aligned with the basic oscillation if they were to effect a dynamical transition. Thus, 4–5 gamma peaks could fit within a single theta positive oscillation, enhancing the underlying oscillation and producing the complex waveform consistent with the experimentally observed beta-gamma/theta nesting.

The second study was a theoretical and computational study of a modified Wilson-Cowan model applied to a 2D lattice of diffusively interacting nodes driven by noise (Di Santo et al., 2008). The lattice was rewired to allow a small fraction of long-range connections. It included a model of synaptic plasticity and some heterogeneity via modified node activity levels. A critical point in the cortical dynamics, indicative of a phase transition in the synchronization, was observed by tuning the baseline level of a synaptic resource parameter. Such a global transition is in contrast to the more localized up-down state transitions found herein. That study also predicts avalanches, as observed experimentally. The synchronization transition was observed in an extended system of size 64²-512² oscillators. Another relevant study is of a system of 80 identical NMs, each connected to 20 neighbors connected in a Watts–Strogatz network, with a specified interaction strength (Kazemi and Jamali, 2022). The system exhibits a synchronization transition as coupling strength is increased, and the transition point shifts with noise stimulus strength.

In each of those studies, the system has a large number of interacting entities, and the transition is driven by interaction strength, while the present study was of a local cluster of six oscillators coupled by the fixed experimentally specified interaction weights and driven by wave-like stimuli. The small size is far removed from the thermodynamic limit of large systems where phase transitions and critical phenomena are usually observed. The transitions observed herein are between up-down states of high and lower amplitude. It would be worth searching for suitable measures of synchronization and their transitions in small groups of NM oscillators and in the pFC cluster. Small groups of oscillators are known to synchronize, thus encouraging such a search. For example, the synchronization of two pendulum clocks has been known for centuries [reported by Huygens in February 1665 (Strogatz, 2003, p. 106)], suggesting that synchronization may be possible in small clusters of oscillators. Whether transitions occur is less certain. The roles of localized up-down states, larger scale synchronization transitions, and the interplay between system size, variability of interactions, and stimuli, all warrant further investigation. The global model of di Santo et al., drawing on the statistical physics of critical phenomena, is more sophisticated than the localized clusters of simple oscillators and so provides a larger perspective on the present study.

The classic WC/JR model has been extended to incorporate the distribution of synaptic weights, leading to a sharper voltage-rate transformation function (Equation 6; Supplementary Equation S5). Heterogeneity has been included via experimentally available linkage weights. Together, these induce interesting dynamical behavior.

The heart of the WC/JR model is a simple, critically damped oscillator wholly driven by external inputs, so it is amenable to intuitive physics to guide its further development and interpretation. This makes clear that brief stimuli, such as gamma bursts, need to be phase-aligned with the underlying oscillations for the stimulus to deliver sufficient energy to drive the oscillations to a new dynamical state. Such phase alignment is reminiscent of the observed nesting of beta/gamma-theta oscillations (Giresch and Plenz, 2008).

A key limitation of NM models is that a number of scales are unclear (cf. Supplementary material 2.7): the spatial extent of a neural mass, the relationship between linkage weights and the synaptic weights, and the relationship of voltage output to measured LFP. The model developed herein can benefit from the inclusion of a model of synaptic plasticity (Tsodyks and Markham, 1997; Di Santo et al., 2008), e.g., via a generalization of the feed-forward and feedback constants A, B, and C (cf. Equations 1–5, Supplementary material 2.2). A biologically plausible way of incorporating synaptic plasticity in the WC/JR model, that avoids arbitrary parameterization, may yield informative dynamical responses. It is unclear if small clusters of NMs can produce synchronization transitions (cf. Section 4.3); more work would be required to establish the cluster size that can produce such.

This study highlights several areas in and near pFC that have interesting roles in local dynamics. Area A32V, the major network in-hub in the marmoset cortex, has a distinctive response to stimuli that warrants further investigation. A11, a connector node and marginal in-hub, is a target of numerous sensory inputs. A46D was sensitive to reset stimuli, indicating a possible role in modulating pFC function; it receives medium-strength inputs from the auditory cortex. The availability of high-density electrode arrays that may be able to probe LFP at the single area level offers an opportunity to study these areas in detail. The new structural data for marmosets at a finer spatial resolution (Watakabe et al., 2023) provides an opportunity to extend the present study.