Brain signals such as Electroencephalogram (EEG), Steady State Visually Evoked Potentials (SSVEP) and Local Field Potentials (LFP) are commonly used in research to understand and identify their correlation with the functional and behavioral states of the brain. Frequency domain analysis is an established, and by far the most popular way of analyzing these brain signals. Broadly, brain signal frequencies are clustered into five bands viz. Delta [δ: 0.5–4 Hz], Theta [θ: 4–8 Hz], Alpha [α: 8–12 Hz], Beta [β: 12–30 Hz] and Gamma [γ: >30 Hz] (1); in many studies, each band is further subdivided into sub-bands (2). Psychophysical studies associate these frequency band oscillations, often termed as “rhythms,” with different brain states in healthy adults. The δ frequency band is usually linked to deep sleep (3); θ oscillations in the neocortex are reported in memory and learning processes as well as in synaptic plasticity (4); α oscillations are most prominent during awake and resting state with eyes closed (1); β oscillations are associated with a high arousal state of the brain for example during sensorimotor activity (5); γ oscillations are associated with attention, consciousness, memory and perception (3). At the same time, alterations in brain rhythms are reliable neuromarkers in several neurological conditions. For example, a reduction in the β band power is a consistent observation in Schizophrenia patients (6); α rhythm slowing is a definitive marker in Alzheimer's (7) and chronic pain patients (8). The differences in the neuronal microcircuits influenced by their underlying synaptic connectivities and attributes are thought to underpin brain rhythmic changes (3); for example decrease in the white matter volume as well as in the overall GABA-ergic inhibition is reported from post-mortem studies of Schizophrenia patients (9). Although much is now known about the correlations between various brain oscillatory states and inter-region connections in the brain, it is hard to understand and identify exact pathway aberrations in disease condition, thereby making it difficult to provide effective treatment. In this regard, in silico models are used by many research labs to study and suggest effective treatment for brain disorders (10, 11). We have been using in silico models as a means to understand better the neural dynamics that underpin various brain oscillatory states. Here, we have used an enhanced version of an existing in silico model to understand the effects of synaptic connectivity changes on the oscillatory dynamics within the θ, α, and β frequency bands.

Recently, we have demonstrated the phenomenon of phase entrainment with periodic flicker visual input in an in silico model (12). The model consisted of three neural structures that lie in the brain visual pathway, viz. the lateral geniculate nucleus (LGN) of the thalamus that receive direct input from the retinal spiking neurons; Layer 4 (L4) of the primary visual cortex that are direct recipient of the visual information from the LGN; Layer 6 (L6) of the primary visual cortex that receive feed-forward input projections from the L4, and provide feedback to the LGN. This layout of the model is an abstraction of the thalamocortical loop in the brain that is essential for generating brain signals for example sleep rhythms (13). The intra- and inter-layer synaptic connectivity parameters for L4 and L6 in the model were based on the anatomical data from the cat primary visual cortex (14). However, the L6 also receives feed-forward projections from the LGN (15), which was not included in the model. In this work, we have presented an enhanced version of the previous model. First, we have included the feed-forward pathway from the LGN to the L6. Second, we have considered all synapses formed between any two populations in L4 and L6 irrespective of the layer where they were formed. For example, the excitatory projection neurons that have their perikarya in L4 also have dendritic projections in layer 5 (L5); similarly, the excitatory projection neurons that have their perikarya in L6 also project their axons to L5 and make pre-synaptic contact with L4 dendrites. This is unlike in our above-mentioned previous work (12), where we considered only those synapses which were made within the L4 and L6. Our modified approach is more biologically realistic and has resulted in an alteration of our base synaptic parameter values compared to that in the previous work. Third, we have addressed missing values in the anatomical data, the details of which are mentioned in section 2.

Our objective in this work is to study the synaptic connectivity changes that underlie alterations of the model signal power in the frequency spectrum for healthy, as well as neurological conditions. We have demonstrated the effect of “lesioning” (disconnecting) synaptic connections in specific pathways of the in silico model that shift the dominant oscillatory frequencies of the neuron populations. We have identified specific pathways in the intact visual thalamocortical loop that cause transitions from the base oscillatory frequency of the neural populations to a state with maximum power within the θ, α, and β bands. The identified synaptic pathways that underlie these state transitions agree with existing anatomical and physiological studies; we have discussed the details in section 4. To test the in silico model for simulating neurological conditions, we have demonstrated the decrease in β band power reported in Schizophrenia patients compared to controls (6, 9, 16). The critical factors to simulate the condition of Schizophrenia in silico are: an overall decrease in the excitatory and inhibitory membrane conductance; reduction in the feedforward excitatory projection signals from the LGN to L4. In summary, our results demonstrate the usability of a simple in silico model of the brain thalamocortical loop, informed by anatomical and physiological synaptic layout and connectivities, to understand the neural pathways that underpin brain oscillatory changes, as well as to validate and inform psychophysical studies.

The layout of this article is as follows: In section 2, we detail the biological basis of the model and the anatomical sources that inform the parameterization of the model synaptic pathways. The simulation methods are also specified. The simulation results and observations are mentioned in section 3. In section 4, we make a detailed discussion of our observations in context to findings from anatomical, physiological and psychophysical research. Concluding remarks are mentioned in section 5.

The thalamocortical model is shown in Figure 1. In section 2.1, we present the biological basis of the model synaptic layout and the data sources for the synaptic connectivity parameters used in the model. In sections 2.3 and 2.4, we present the parameterization and simulation methods, respectively, of the model.

The importance of understanding the synaptic connectivity correlates of brain oscillatory changes cannot be more emphasized than in neurological disorders (25, 26). As brain rhythms are impacted during neurological conditions, they become an important neuromarker for the first line of clinical investigation (27). For example, slowing of the frequency spectra recorded in EEG is a definitive neuromarker of Alzheimer's disease, a condition inflicted by aggressive loss in brain synaptic connectivity (2, 7); dysfunctional beta oscillations due to anomalies in the networks of GABA interneurons and in cortico-cortical connections are associated with cognitive deficits in Schizophrenia (9, 28, 29). Therefore, it becomes critical to understand the exact synaptic pathway correlates of such neuromarkers, and possible therapeutic solutions. One clinically successful instance of therapeutic solution for the EEG neuromarker (unwanted 10 Hz oscillation) of Parkinson's disease is Deep Brain Stimulation (DBS) (30). However, clinical investigation of possible therapeutic treatments is hard to conduct in neurological patients; this is especially problematic in instances of cognitive impairment, where patients may feel distress, which in turn may aggravate their condition. Animal models of neurological conditions have provided a lot of clues, but human and primate brains react differently than animals, and findings may be misleading. There are other organoid solutions being suggested recently (31). One approach that is relatively cheaper and have low demands on resources are in silico models that are brain inspired. Thus, these models can be used to represent neural circuitry governing brain rhythms and deepen our understanding of healthy as well as sick brain states. They also come as an important aid in exploring new, as well as enhancing existing, therapeutic treatments for neurological conditions (10, 11).

We have been using in silico models of the thalamocortical loop, informed by anatomical and physiological studies of the brain visual pathway, to simulate brain rhythms within the alpha and theta band frequencies. Such frequency domain analysis of brain signals such as EEG, SSVEP and LFP are by far the most popular means to visualize and correlate specific frequency sub-bands [viz. delta (δ), theta (θ), alpha (α), beta (β), and gamma (γ)] to brain states in both health and neurological conditions. In a recent work, we have validated a model with psychophysical studies that demonstrated phase entrainment of brain signal during SSVEP. In this work, we have presented an enhanced version of this model by adding more biologically plausible synaptic layout and connectivity parameters. Our in silico model consists of three modules viz. the lateral geniculate nucleus (LGN), that is the visual thalamus receiving direct input from the retina; Layers 4 and 6 of the primary visual cortex (V1), both of which receive feedforward input from the LGN; the L6 provide feedback to the LGN. The model synaptic layout is as in Sherman (15). The synaptic parameters for intra- and inter-layer connectivity in L4 and L6 are computed from Binzegger et al. (14). The intra-layer connectivity in LGN and feedback from V1 layers are as in our previous works (2). For missing data in our sources, informed assumptions are made. The voltage time series outputs of the model neuronal populations are converted to the frequency domain to identify the dominant frequency of oscillation. The outputs of the excitatory projection neurons of each module, viz. the Thalamocortical Relay (TCR) cells of the LGN, and Pyramid cells of L4 (Py4) and L6 (Py6) are observed for demonstrating brain state transitions. Where we have validated our results with clinical data, the output of the L4 is considered. This is in agreement with the observations made in Maier et al. (19), where the current source density (CSD) of LFP recordings from V1 following a visual stimulus had origins in L4 and was attributed to the thalamocortical afferents.

One approach for investigating clinical conditions with in silico model is to set a “control” state that will simulate the specific brain state (which the model is intended to investigate) in healthy adults. This control state acts as a reference point from where model parameters and attributes that produce aberrations corresponding to the disease states can be identified. We have adopted a similar approach in this work. At first, we recorded the base state oscillations of the model that emerged when the synaptic connectivity parameters in the model were set to the base values mentioned in Table 3. Next, starting from the base state of the model, we simulated transitioning into the α, β, and θ states, and identified the underpinning synaptic parameter pathways. Note that for each state, we aimed at all the excitatory projecting populations, viz. the TCR of LGN, Py4 of L4 and Py6 of L6 to have a peak power within the respective bands. This is because, to the best of our knowledge, there is no coherent data in literature that suggest the possible oscillatory combinatorics in LGN and V1 and their associated brain states. Furthermore, literature suggest that same areas of brain may employ multiple mechanisms to generate multiple frequency bands based on the brain state and function (32). We too have observed an exploding parameter and oscillatory state combinatorics in the model. We have chosen to demonstrate the synaptic connectivities that can produce the same peak oscillatory frequency in the TCR, Py4 and Py6 populations simultaneously.

In the base state of the model, the peak frequency of oscillation for TCR and Py6 is in the α band, and for Py4 is in the θ band. For long, the thalamus has been attributed as the primary pacemaker of α oscillations (33, 34). However, LFP recordings from the V1 of macaque monkeys and dogs suggest α oscillation origins in the deep cortical layers like Layer 5 (L5) (35, 36) and L6 (37). Another LFP based study of the occipital lobe of macaque monkeys suggest the superficial cortical layers 2/3 (38) as the origin of α rhythms. Bollimunta et al. (37) and Mila et al. (38) also discuss the possibility of separable thalamic and cortical α pacemakers which become differentially active and coupled under different behavioral conditions. Interestingly, two mesoscopic computational models where the α oscillations can be generated with or without the thalamic drive is reviewed in Sigala et al. (39). To understand the origin of the α rhythm in our base state network, we lesioned the thalamus from the V1 layers by disconnecting all synaptic feedforward and feedback pathways from and to the LGN, respectively. All other connections were left intact. This lesioning did not affect the α oscillations in Py6; rather, the TCR α peak moved from the α high (α-H) to the α low (α-L) sub-band due to lesioning the Py6 feedback. Thus, our model supports the discussions in Bollimunta et al. (37), Mila et al. (38) where both LGN and the infragranular L6 in V1 demonstrate α rhythm generation independently. In future works, we have plans to integrate the infragranular L5 and supragranular L2/3 in our model, which will allow us to test this aspect more thoroughly.

In addition to the above, lesioning of the LGN also did not affect the Py4 peak frequency within the θ band. To understand the origins of this θ oscillations, we lesioned the L4 and L6. The Py4 peak frequency shifts to the α band, but the Py6 peak frequency moves from α-L to α-H. Investigating further, we identified that the inhibitory feedback from Py6 to Py4 is the primary cause of the θ dominance in Py4. In turn, the feedforward excitation of L4 on L6 slows down the α band power in Py6.

Next, from the intact base network state, we aimed at identifying synaptic connectivity parameters that could cause brain state transitions to θ, α, and β bands for TCR, Py4 and Py6 simultaneously. Our findings are as follows:

This concludes our study of the control state conditions of the model in the θ, α, and β band oscillatory state.

As a case study, we demonstrate the known alterations in θ and β band frequencies in the condition of Schizophrenia. It is reported by Adams et al. (6) that resting state EEG of Schizophrenia patients show a decrease in β band power and increase in θ band power compared to controls. (They have also mentioned alterations in the γ band, but we do not consider the γ band oscillations in this work). At the same time, the authors make a computational study to validate their psychophysical results and propose possible synaptic correlates of the frequency band alterations in Schizophrenia conditions. Overall, they identified several combinations of excitatory-inhibitory imbalances that simulated the β and θ band power alterations in their model. Similar combination of psychophysical and animal model-based studies also indicate the excitation-inhibition imbalance as the neuronal correlates of the β band alterations observed in Schizophrenia patients compared to controls (9). We observed that two critical changes in our model induced a significant decrease and increase (with P < 0.001) in β and θ band power, respectively: (a) decrease in the AMPA and GABA_A synaptic conductance parameter values; (b) decrease in the excitatory synaptic connectivities from LGN to L4. The observations in (a) conform to those in the above-mentioned researches (6, 9). The observations in (b) are particularly interesting, where visual pathway impairment is identified by many researches as a biomarker of Schizophrenia (44). Along these lines, an SSVEP study by Giri et al. (45) demonstrated significant decrease (P < 0.001) within the β band power in SSVEP outputs corresponding to selected visual flicker stimuli. Tested with our model for simulated conditions of Schizophrenia, the results agree with these observations.

Readers may note that the condition of Schizophrenia in our work is simulated with the β band oscillatory model state as the control condition. However, in Adams et al. (6), the EEG studies were made when all controls and patients were in the awake resting state and showing an α peak oscillation. While we could generate a reduced β band power when the model is in the α oscillatory control state, we did not see an increase in the θ band power, thus unable to fully validate the experimental results in the model. The main criticism of the model presented here is that the synaptic as well as the structural layout of the V1 module is simplified. Also, synaptic connectivity parameters in the model are informed by the cat visual cortex, but the Schizophrenia symptoms that are being validated are observed in humans; this is in spite of the primate visual system having evolved with significant differences compared to cats. As a future work, we will enhance the V1 model by adding the superficial L2/3 and infragranular L5 and parameterizing based on primate anatomical and physiological data (46). Furthermore, we aim to include γ rhythm observations in the model, as these are also identified as possible neuromarkers in the EEG of Schizophrenia condition. Also, we plan to use cross-frequency coupling (CFC) measures to analyse our model dynamics, as CFC allows a better understanding of thalamus and neocortex information processing (47) and abnormal behavior in neurological diseases (48).

The in silico model presented here belongs to the neural mass genre of computational models (NMM) that can simulate complex biological phenomena with fewer state variables and parameters. This make NMMs computationally less expensive compared to spiking neural networks (SNN), although at the cost of single-neuron-level dynamics unlike in the latter (49). In addition, the sigmoid function used in NMMs to represent the neuronal activity for a population is not based on a detailed biophysical description of spiking neuron (50). Subsequently Byrne et al. (50) suggested a change in NMM tactics and proposed a next generation NMM. Recently, Huang and Lin (51) proposed a current density based model of NMM. Interestingly, Ruffini et al. (52) proposed a laminar based NMM for transcranial current stimulation (tCS). tCS is emerging as a promising alternative treatment in drug resistant Schizophrenia patients (53) and hence modeling this use case gives opportunity to personalize the treatment for individual patients. As future work, we see potential to extend the Schizophrenia case-study presented here to better understand the effects of tCS in the condition.

In conclusion, we have tested a simple in silico neuron population network of the thalamocortical loop in the brain visual pathway. The model could generate brain rhythms within the α, θ and β bands and correlate with specific synaptic connectivities and pathways. As a case study, we have simulated the oscillatory power alterations as observed in experimental studies on Schizophrenia patients. Our results show a significant (P < 0.001) decrease in beta band power and a simultaneous increase in the theta band power, similar to that observed in Schizophrenia.

The main contribution of this work is a thorough analysis of the synaptic connectivity parameter combinatorics that generates the α, θ, and β oscillatory states in the in silico model; for brevity, we have presented only selected combinations in this work. Findings from anatomical and physiological studies suggest that all three oscillatory states are observed in both LGN and the V1. It is however hard to identify the exact synaptic connectivity combinations that underpin the continuum of changing brain states in vivo. Our findings demonstrate the potential for designing and using simple brain inspired neural networks to simulate both control and neurological conditions. This in turn has the potential to be used for investigating individual therapeutic requirements to reverse or suppress those conditions. As ongoing work, we are developing this model to allow deeper insight into the synaptic connectivities and attributes that underpin the conditions of Schizophrenia.

The original contributions presented in the study are included in the article. Further inquiries can be directed to the corresponding author/s.

SS: methodology, software, investigation, writing—original draft, and writing—review and editing. BS: conceptualization, methodology, writing—original draft, writing—review and editing, and funding acquisition and supervision. Both authors contributed to the article and approved the submitted version.

This research was supported by the Birla Institute of Technology and Science (BITS) Pilani Goa Campus Grants BPGC/RIG/2018-19 and GOA/ACG/2019-20/Oct/02 to BS. SS was supported by the BITS Pilani Goa Campus Institute Fellowship awarded toward her Doctoral Research.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.