In silico hippocampal modeling for multi-target pharmacotherapy in schizophrenia.

Treatment of schizophrenia has had limited success in treating core cognitive symptoms. The evidence of multi-gene involvement suggests that multi-target therapy may be needed. Meanwhile, the complexity of schizophrenia pathophysiology and psychopathology, coupled with the species-specificity of much of the symptomatology, places limits on analysis via animal models, in vitro assays, and patient assessment. Multiscale computer modeling complements these traditional modes of study. Using a hippocampal CA3 computer model with 1200 neurons, we examined the effects of alterations in NMDAR, HCN (Ih current), and GABAAR on information flow (measured with normalized transfer entropy), and in gamma activity in local field potential (LFP). We found that altering NMDARs, GABAAR, Ih, individually or in combination, modified information flow in an inverted-U shape manner, with information flow reduced at low and high levels of these parameters. Theta-gamma phase-amplitude coupling also had an inverted-U shape relationship with NMDAR augmentation. The strong information flow was associated with an intermediate level of synchrony, seen as an intermediate level of gamma activity in the LFP, and an intermediate level of pyramidal cell excitability. Our results are consistent with the idea that overly low or high gamma power is associated with pathological information flow and information processing. These data suggest the need for careful titration of schizophrenia pharmacotherapy to avoid extremes that alter information flow in different ways. These results also identify gamma power as a potential biomarker for monitoring pathology and multi-target pharmacotherapy.

Introduction

Schizophrenia is a chronic disease with a lifetime prevalence of around 4/1000[1], which usually produces life-long disability[2]. Cognitive impairment and information processing deficits are chief causes of disability[3-5]. The most affected cognitive domains are processing speed, working memory, episodic memory, and verbal learning and memory[6-8]. Current antipsychotic medications have limited impact on cognitive symptoms and information processing deficits[9,10]. Therefore, there are significant gaps in therapy and patients’ clinical care[11].

Recent research has emphasized the role of glutamatergic transmission as an extension of the dopaminergic hypothesis for schizophrenia pathophysiology, especially to capture cognitive impairment associated with schizophrenia (CIAS) and information processing deficits[12]. The role of glutamatergic transmission has been supported by the psychotomimetic effect of N-methyl-D-aspartate (NMDA) receptor (NMDAR) antagonists like phencyclidine (PCP) and ketamine in healthy volunteers[13,14]. Ketamine also worsened cognitive symptoms in schizophrenia patients[15].

Oscillations are abnormal in schizophrenia patients, both in the low (delta and theta) and high (gamma) frequency ranges[16,17]. For example, a positive correlation was found between power of theta oscillations in the temporal lobe and positive symptoms using resting state MEG[18]. In a 40 Hz auditory steady-state response (ASSR) paradigm, theta power was found to be increased, and this was associated with impaired verbal memory[19]. Cortico-hippocampal theta oscillations play multiple roles in memory encoding and retrieval (reviewed in ref. [20]), and mediate sequencing of events through theta/gamma coupling[21,22].

Gamma oscillations are also involved in sensory integration[23], enabling local computations within cortical microcircuits[24] by providing binding among neurons belonging to a particular neuronal ensemble[25,26], and routing information across brain regions[27]. Multiple studies revealed reduction in induced and evoked gamma power[28-31] in schizophrenia patients. By contrast, presymptomatic, clinically high-risk individuals demonstrated an increase in resting gamma power[32]. Ketamine increased resting gamma power in healthy volunteers[33]. Interestingly, acute ketamine increased resting and evoked gamma power in rodent hippocampal CA3[34], while chronic ketamine reduced gamma power[35], similar to what has been reported in patients with chronic schizophrenia[32].

Combining evidence from gamma studies and ketamine studies suggests a role for NMDAR dysfunction in gamma abnormalities and in CIAS. Additional evidence comes from the schizophrenia genome-wide association study (GWAS) published in 2014, which identified GRIN2A (glutamate ionotropic NMDA-type receptor subunit 2A) on chromosome location 16p13 as being close to one of the 108 gene loci identified[36]. Clinically, agonists for NMDARs that have been studied for the treatment of CIAS include glycine and D-serine, both of which bind to an allosteric site on the GluN1 subunit of NMDAR and are obligatory for NMDAR activation by glutamate[37]. D-cyclo-serine, a partial agonist at the glycine site[38], has also shown efficacy in some studies. However, these studies yielded mixed results; a recent meta-analysis found no overall significant effect on cognition, although young patients aged 30–39 years old showed some benefit[39].

Another molecular aspect that we consider in this paper is the involvement of the HCN1 gene (hyperpolarization-activated cyclic nucleotide-gated channel type 1) that codes for Ih channels (hyperpolarization-activated current, also known as anomalous rectifier, If, Iq)[40-42]. Along with NMDAR, this channel plays a role in the generation and modulation of neuronal oscillations[43,44]. The HCN1 gene on chromosome location 5p21 is also close to one of the GWAS loci associated with schizophrenia[36]. The role of Ih in oscillations, and evidence from the GWAS study, suggest HCN1 product as another potential therapeutic target for CIAS[45].

Manipulating the GABAergic system provides another potential treatment target[46]. The GABAergic system shapes synchronized neuronal activity during oscillations[47-51], and is impaired in schizophrenia. Postmortem studies have revealed low inhibitory interneuron glutamic acid decarboxylase enzyme (GAD67), low GABA transporter, and low pyramidal cell GABAA receptor (GABAAR) α1 subunit mRNA transcripts in the frontal lobe of patients[52]. Similar findings were also demonstrated in the hippocampus, where a reduction in numbers of somatostatin-positive and parvalbumin-positive interneurons was also found[53]. The GABAergic deficit hypothesis for schizophrenia is further supported by reduced postmortem immunoreactivity of GAD65/67 in interneuronal neuropil in the hippocampus[54]. Reduced GABA tone has been suggested to mediate hippocampal hyperactivity in these patients[55-59].

Schizophrenia pathology spans the multiple hippocampal subfields, as well as affecting adjacent para-hippocampal areas. We focus here on area CA3, an area that shows increased activity in schizophrenia but is relatively preserved in comparison to CA1. CA3 hyperactivity is seen in schizophrenia patients[56], and in animal models of psychosis[60,61]. CA3 projects strongly to CA1, which also shows evidence of hyperactivity[62,63], presumably as a consequence of these projections. CA1 shows shrinkage, even early in the illness[64-66], suggesting early and prominent cell loss in CA1. A recent postmortem subfield transcriptome analysis of schizophrenia suggested that immune-mediated changes in CA3 might be one cause of CA3 excitability, which could subsequently produce synaptic scaling after cell loss in CA1[67], thereby further augmenting CA1 neuron loss[68]. Another reason to focus on CA3 is its possible role in pattern completion[69-71], a component of cognitive processing. In this context, a hyperactive CA3 might alter pattern completion through activity spill-over[59,72,73].

As described above, schizophrenia pathophysiology is extremely complex, with abnormalities at multiple scales, from genes, second messenger cascades, and cells, up to local networks and inter-areal communication. Given this complexity, it is reasonable to consider that multi-target pharmacotherapeutic approaches will be useful[74]. Targets in multi-target pharmacotherapy interact in highly non-linear ways, making it impossible to intuitively predict the effects[75-77]. We, therefore, use simulations to study these interactions. In this study, we investigated how alterations in the three targets mentioned above—NMDAR, Ih, GABAAR—will affect (1) oscillations and (2) information flow in a biophysically-realistic computer model of the CA3 region of the hippocampus. This builds on our prior results which identified the effects of blocking NMDARs on oriens-lacunosum moleculare (OLM) interneurons on gamma oscillations and information flow[78].

Results

Baseline network activity produced theta with nested gamma oscillations

More than 2000 simulations were run using NEURON 7.4[79,80]. These included 175 simulations for the NMDAR augmentation (5 input random seeds × 5 connectivity random seeds × 7 NMDAR scalings), and 1575 simulations for the tri-scaling interactions of NMDAR, Ih, and GABAAR (3 input random seeds × 3 connectivity random seeds × 7 NMDAR scalings × 5 Ih simulations × 5 GABAAR simulations), with additional preliminary simulations to tune the model. Multiple versions of individual simulations with different randomizations were run to ensure that results were reliable and not due to specific network architecture or background input. An individual simulation of 7 s model time required about 5 min computer time using an integration timestep of 0.1 ms on a Linux system with eight 2.67-GHz Intel Xeon quad-core CPUs. The full study utilized ~160 core-hours of high-performance computing time. Simulations were also run with finer timestep to assess numerical stability. We investigated the correlations between ion channel manipulations, scaling , , and , with information flow and gamma oscillation strength. We present the baseline simulations first, then NMDAR augmentation simulations, followed by the results from the simultaneous manipulation of NMDAR, Ih, and GABAAR.

In baseline simulations (Fig. 1), firing of the three neuronal populations was synchronized at two primary frequencies: gamma (~25–50 Hz) and theta (~6–8 Hz). Synchronized activity resulted from the following sequence of events: pyramidal neuronal firing triggered firing of PV basket interneurons, which in turn turned off the pyramidal neurons until inhibition wore off (pyramidal interneuron network gamma, PING oscillations) (Fig. 1c). Pyramidal neurons also triggered the firing of OLM interneurons, which provided further feedback inhibition to the pyramidal neurons, but at a slower rhythm (theta) due to OLM interneurons’ longer-lasting inhibition from longer GABAAR time constants. Synchronization at the theta frequency was also mediated by inhibitory input from the medial septum to both OLM and PV interneurons, but was not dependent on it. While pyramidal neurons were under the inhibitory influence of the OLM interneurons, the firing and reciprocal inhibition between PV basket interneurons generated interneuron network gamma (ING) oscillations (Fig. 1d). The synchronized firing of the pyramidal neurons, PV basket, and OLM interneurons was reflected in the local field potential (LFP) (Fig. 1b). The LFP showed both gamma and theta frequencies. The power spectral density (PSD) (Fig. 1e) showed the peaks of gamma (~35 Hz) and theta (~7 Hz) oscillations in this representative simulation.

Fig. 1

Baseline simulation generated theta (6–8 Hz) with nested gamma (30–40 Hz) rhythm.

a Raster plot where each row represents a neuron and each dot represents an action potential (x-axis: time; y-axis: neuron identity; OLM interneurons in blue, PV basket interneurons in green, pyramidal neurons in red). b Simulated local field potential (LFP). c, d Simplified network diagrams showing mechanisms generating PING and ING oscillations. Arrows indicate activation while bar-headed lines represent inhibition. e Power spectral density (PSD) of LFP showed theta (~7 Hz) and gamma (30–40 Hz) peaks.

Augmenting OLM NMDARs modulated circuit oscillations and information flow

Our first step was studying the effects of augmenting NMDAR of OLM interneurons on gamma oscillations and information flow across the model (Fig. 2). OLM interneurons were the location where NMDAR antagonism produced theta and gamma oscillatory changes similar to what was seen with ketamine in hippocampal CA3[34,78]. Taking the simulations in Fig. 1 as the baseline, we then scaled of OLM interneurons from 1.5× to 30× the control simulations. This scaling could reflect a conductance increase due to phosphorylation or insertion of different isoforms with or without an increase in the number of NMDAR channels. Scaling of OLM interneurons reduced gamma power ~10-fold until it disappeared at above 30× control (Fig. 2b, left y-axis).

Fig. 2

Augmenting on OLM interneurons.

a Information transfer (nTE) and modulation index (MI) showed an inverted-U relationship with scaling OLM up (mean ± SEM). b Increasing OLM increased activity of OLM interneurons (blue), reduced pyramidal (red) and PV (green) neuronal firing, and reduced gamma power (mean ± SEM). c Inverted-U relationship between gamma power and information flow shown with different levels of OLM (mean ± SEM). Note that highest is now at the left (values in shaded rectangles). d Raster plots and LFP at different scaling relative to control.

Scaling OLM upward increased drive onto OLM interneurons and therefore increased the firing rate of these inhibitory cells (Fig. 2b, blue line), which reduced the firing of the other two cell types that received OLM projections (PV, green; PYR, red). Information transfer (nTE) across the excitatory (PYR) population (Fig. 2a black) peaked at an intermediate level of PYR and PV firing, and showed far lower values with firing rate extremes, low or high. Cross-frequency coupling between theta phase and gamma amplitude, measured by modulation index (MI), also showed the inverted-U relationship with NMDAR scaling, correlating with information flow (Fig. 2a red). We found a similar inverted-U that related information flow to gamma power, a measure of synchronized firing (Fig. 2c). Note that the values of are reversed in Fig. 2c relative to 2a and 2b, with high on the left associated with low PYR firing and low gamma. There were two causes for low gamma with low PYR rate: (1) reduced number of spikes, which also reduced the strength of theta; (2) reduction in spiking coherence which was no longer as sharply shaped by PING and ING (Fig. 2d, scale = 30).

Modulating information flow with multiple channel manipulations

Consideration of multi-target pharmacotherapy requires determining how actions at different targets combine. We, therefore, looked at alterations in Ih and GABAAR, both of which are considered possible factors in schizophrenia pathophysiology[36,53], along with . Having previously demonstrated strong effects on gamma of HCN channels of pyramidal and PV neurons[43], we focused here on Ih at these two locations. The inverted-U configuration of the nTE peak for augmentation was seen at all Ih scalings except for 10× (Fig. 3). An inverted-U pattern could also be seen for GABAAR scaling with 10×, again at all but the highest Ih value. Similarly, the inverted-U peak can be detected for Ih scaling around the point marked H1, with a less well-defined, broader surface peaking around the point marked by H2. Thus, we have identified parameter ranges that identify a set of peaks which are points of strong nTE in a 4-D space based on this 3D parameter space.

Fig. 3

Multi-target manipulations modulate information flow.

Ih scaled in panels (a–e), with and on x- and y-axis in each case. Values are interpolated from simulation results (small circles). L1 and L2 are points in the tri-scaling parameter space with low information flow, while H1 and H2 are points in the tri-scaling parameter space with high information flow.

Firing patterns from regions of low (L1–2; Fig. 4) and high (H1–2; Fig. 4) nTE were easily distinguishable. Two patterns of neuronal firing were seen which produced low nTE: (L1) low firing; (L2) over-synchrony. L1 exhibited low oscillatory power in addition to low firing. This was because of the combined inhibitory effect of high NMDAR scaling on OLM interneurons (30×) and high GABAAR scaling on pyramidal neurons (10×), leading to inhibition of pyramidal neurons and PV basket interneurons. In L2, highly synchronized phase-locked pyramidal neuron firing was driven by strong PV basket feedback inhibition, permitting little of the variation required to transfer information. This highly synchronized dynamics resulted from high Ih scaling (10×) on pyramidal neurons and PV interneurons, increasing their excitability and their synchronous activity, similar to what we reported before[43]. The different dynamics between L1 and L2 were reflected in their gamma power (bottom panel of Fig. 4).

Fig. 4

Intermediate excitability with intermediate synchrony allowed for high information transfer across tri-scaling manipulations.

Raster plots based on parameters labeled L1, L2 (with reduced information flow), and H1, H2 (with increased information flow) on Fig. 3 for an example simulation. L1 showed low firing; L2 showed high synchrony. H1 and H2 showed similar firing patterns, where both excitability and synchrony were midway in comparison to simulations in L1 and L2. In bottom panel, gamma power reflected the different raster dynamics, with intermediate gamma power for H1 and H2 simulations in comparison to gamma power of simulations L1 and L2. nTE showed an inverted-U pattern with increasing gamma power.

Examination of rasters from points of high information flow-through (H1–H2 in Fig. 3) showed moderate levels of activity (H1–H2 in Fig. 4). Firing here was intermediate between the low firing of L1 and the strong phase-locking of L2, consistent with the inverted-U pattern. This intermediate synchronization allowed the pyramidal neurons to be excitable enough for the timing of their firing to be appropriately biased by driving input to allow information flow. Gamma power of H1 and H2 in the bottom panel of Fig. 4 was also intermediate between gamma power of L1 and L2. This pattern generalized across all manipulations, where maximum nTE was found at a mid-gamma-power range. As shown in Fig. 5, gamma power from the tri-scaling simulations was in the range of 0.01–10.8 mV2/Hz. However, the simulations with nTE in the top 10th percentile had gamma power at an intermediate range of 0.02–2.12 mV2/Hz.

Fig. 5

High nTE values were found at intermediate gamma power, suggesting that extremes of gamma power reduced information flow.

Top panel shows gamma power versus nTE in simulations from tri-interaction pharmacological manipulations. Each dot represents a simulation. The red set represents manipulations with nTE values in the top 10th percentile, which lied in an intermediate gamma-power range from 0.02 to 2.12 mV2/Hz. Bottom panel shows a histogram of the distribution of the data points in the top panel.

The interaction between excitability and synchrony could be seen in the inverted-U shaped relationship of information flow with the firing rate of pyramidal neurons and PV interneurons (Fig. 6). At low-firing rate of pyramidal neurons, reflecting low excitability, there was reduced information flow (Fig. 6a, left). As firing rate increased, information flow increased. However, after a certain degree of excitability, information flow started decreasing. This could be explained by the limiting effect of highly synchronized firing on information flow, as reflected in the population synchrony of pyramidal neurons (red markers in Fig. 6c). Population synchrony was estimated using the coefficient of variation from the interspike intervals across all the neurons in a population[81]. Population synchrony ~1 represents interspike intervals that belong to a Poisson distribution, reflecting random population firing driven by the Poisson driving input. In a synchronized population, the upper range of population synchrony approaches , where N is the number of neurons firing synchronously at each firing cycle. At highly synchronous firing of pyramidal neurons (red markers in Fig. 6c) and PV interneurons (red markers in Fig. 6d), information flow was reduced. The coloring in panels a–d of Fig. 6 is for quintiles of pyramidal neuronal population synchrony, to keep track of the simulations with similar manipulations across firing rate and population synchrony. There was no clear relationship between OLM interneuronal firing pattern and information flow (Supplementary Fig. 1).

Fig. 6

Pyramidal and basket cell population firing rates and synchrony were associated with information flow changes.

Variability of population firing measured with population synchrony. Points in the panels, representing individual simulations, are color-coded according to quintiles measured for pyramidal population synchrony. Quintiles show a general consistency across these four measures.

Discussion

Using a computer model of the CA3 microcircuit, we found an inverted-U pattern which related information flow both to potential therapeutic targets for schizophrenia—NMDAR, GABAAR, Ih, and to gamma activity. These findings suggest an interesting interaction between pyramidal neuronal excitability and synchrony (Fig. 7). At low excitability, pyramidal neurons were below firing threshold and showed low activity. A driving input was less likely to alter spike timing. Therefore, information flow was low. As excitability increased, more pyramidal neurons were closer to the firing threshold, and the driving input could more strongly bias firing, increasing information flow. At the other extreme of excitability, high pyramidal neuronal activity increased activity of basket interneurons, increasing gamma power via both ING and PING mechanisms (Fig. 1). With increased gamma power, pyramidal cell firing was increasingly locked into the oscillation, reducing the ability of inputs to influence firing times and reducing information flow. These findings provide a mechanistic explanation connecting pyramidal neuronal excitability and activity, with population dynamics and oscillations, in the context of potential therapeutic targets in schizophrenia.

Fig. 7

Interaction between pyramidal neuronal excitability and synchrony affected response to driving input, and thus information flow.

Driving inputs arriving at a population of pyramidal neurons (red circles) increases excitability (arrows). At low excitability, driving input is not enough to reach threshold and trigger firing, reducing information flow from driving input. At moderate excitability, pyramidal neurons population is close to firing threshold and so driving input is enough to push cells into firing. At high excitability, pyramidal neurons are pushed back-and-forth between synchronized firing with little driving input influence relative to internal drive, and synchronized inhibition with little input influence due to distance from threshold (RMP: resting membrane potential).

Our study makes the following testable predictions: (1) Augmenting NMDARs on OLM interneurons, for example using photo-switchable NMDAR GluN2A or GluN2B subunits[82] expressed specifically in the OLM interneurons, will show alterations in gamma activity in vivo. (2) The effect of augmenting NMDARs on information flow will depend on the dynamic state of the network reflected by gamma power. Optogenetic stimulation of channelrhodopsin-2 expressed in the PVs at gamma frequency could be used to manipulate gamma power[83]. Information flow could be measured by nTE using LFPs[84]. (3) Manipulations of Ih with photo-activated adenylate cyclase expressed specifically in PVs and PYRs, or by using norepinephrine, will produce our observed changes in gamma and information flow.

Either extreme of information flow will interfere with information processing in the brain. Brain organization is thought to be grossly hierarchical[85]. As a first approximation, we can ignore feedback circuits and consider the processing of information at each level from the prior processing stage[86,87]. We have previously hypothesized, and demonstrated in a model, that proper information processing at each stage requires a balance between information flow-through and information from the area itself, related to local resonance properties[88]. In the present model, where there is relatively low information flow-through, flow-through maximization would be optimal. Decreased flow-through across the CA3 microcircuit would provide decreased information at CA1, the next step in processing, with activity primarily related to CA3 dynamical state rather than to input from dentate gyrus.

Excessive dependence on internal activity has been proposed as a possible pathophysiological mechanism underlying the development of delusions in schizophrenia. Similar to our conclusions, Tamminga et al.[59] suggested that CA3 hyperactivity in schizophrenia patients might indicate that the hypothesized pattern completion function of CA3 is excessive, resulting in faulty pattern completion that is not based on environmental input, producing delusions. In the cognitive domain, decreased information flow would impair performance because information being conveyed would be related more strongly to the dynamical patterns intrinsic to the CA3 microcircuit, rather than being related to the current task or to sensory input.

Coordination between different neural processes (neural coordination), including neuronal firing and oscillations, has been proposed to be important in cognitive coordination required to coordinate two or more frames of reference (e.g., visual and olfactory frames of reference in rodents), whether maintained in separate ensembles simultaneously or provided by alternating ensemble activation at delta or theta frequency[89]. Neural discoordination would then underlie cognitive symptoms (cognitive discoordination) in schizophrenia[16,90,91]. Failure to coordinate multiple frames could result in difficulty separating frames reflecting internal processing (expectations, Bayesian priors, imagining, planning) from frames reflecting external realities (stimulus associated). In our model, moderately-synchronized pyramidal neuronal firing was needed to increase information flow.

Impaired excitation/inhibition balance has been proposed to be one of the pathophysiological mechanisms in schizophrenia. In a combined fMRI and computer modeling study, Starc et al.[92] showed impaired working memory deficits suggestive of cortical disinhibition. Impaired PV functioning in schizophrenia suggested the shift of the excitation/inhibition balance toward more excitable microcircuit as well[93-95]. Our model suggested how more excitable microcircuit can reduce information flow through increased synchrony.

Other groups have studied the relationship between gamma power and information flow outside the context of schizophrenia. Buehlmann et al.[96] showed increased transfer entropy with increased gamma power between two sets of excitatory and inhibitory integrate-and-fire neurons. In addition to having a simpler model, they performed a modulation different than what we used here, the ratio, which may not have permitted them to explore as broad a range of gamma power. Akam et al.[97] described a model where oscillations provided an additional channel for information flow besides the rate of neuronal firing. They studied the flow of information from four input populations of excitatory and inhibitory neurons into one output population. Switching the dynamics of one input population from irregular to oscillatory firing resulted in additional information flow, consistent with our study.

The major limitations of this study are the limitations that are inherent in all modeling studies—we necessarily made choices as to what to include and what to leave out. The things left out include both “known unknowns ... and unknown unknowns.” With respect to the unknown unknowns, we have limitations that are comparable, though more severe, than those of other studies—once you move to an animal disease model (in vivo) or remove tissue from the organism (in vitro) or simply extract parameters (in silico), you are eliminating much of the clinical phenomenology. With respect to the known unknowns, we are progressively adding detail and specifics using the best information available, but we continue to have computing and research limitations that reduce the fidelity of the model. In particular, (1) we omitted interneuron populations other than PV, OLM cells; (2) we omitted dopaminergic and serotonergic receptors, the targets of most current antipsychotic medications; (3) we modeled PV and OLM interneurons as having the same proportions; (4) we modeled PV and OLM as single compartments without 3D details; (5) there is inadequate information regarding the distribution of voltage- and calcium-sensitive ion channels in PYR dendrites.

Our multiscale modeling study suggests that network synchronization and pyramidal neuronal excitability are potential dynamical targets for the treatment of cognitive symptoms and information processing deficits in schizophrenia. We showed how manipulating multiple molecular elements in a multi-target pharmacotherapy approach explains some of the inverted-U shaped phenomena seen in schizophrenia[98]. Therefore, each patient may have to be carefully “tuned” to the middle range (either with DBS or pharmacology). Such an approach, a personalized medicine approach, fits the heterogeneity of the symptoms and genetic risk for the illness.

Multi-target pharmacotherapy is necessarily already used in schizophrenia treatment. Clozapine, often used after the failure of multiple other antipsychotic medications[99], is considered a “dirty” medication that targets a wide array of receptors[100]. A multi-target pharmacotherapy approach is also supported by the evidence that multiple genes and proteins have been identified in the pathophysiology of schizophrenia (e.g., see ref. [101]).

A number of FDA approved pharmacological agents act on the molecular targets we investigated here, and so could be explored for multi-target drug therapy complementing NMDAR augmentation in the treatment of cognitive symptoms in schizophrenia. GABAAR modulators include benzodiazepines, which are not subunit-specific, as well as subunit-specific agents, e.g., zolpidem[102]. But until now, GABAAR modulators are not neuronal subpopulation specific. In regard to Ih, two medications, lamotrigine and gabapentin, used for various neurological and psychiatric disorders, upregulate HCN1 levels[103].

Our study points to a possible explanation of a non-intuitive relationship between gamma power and symptoms of schizophrenia. In schizophrenia and pharmacological models of schizophrenia, gamma power has been found to be decreased[16,104,105] or increased[32,106,107], depending on the study. Our results suggest a way to resolve the paradox: both extremes of gamma-power interfere with information processing. Psychotic symptoms, such as delusions of control and hallucinations, as well as cognitive symptoms, have been conceptualized as being due to “dysconnection” syndromes[108], where the communication between different brain regions is disrupted. We would instead suggest that the complexity and variability of psychotic manifestations might instead be due to re-connection or re-wiring syndromes, where areas of decreased oscillatory strength will be disconnected and areas of increased oscillatory strength hyper-connected but with little variability and responsiveness to outside inputs (stereotypic thoughts and behaviors). This rewiring would be due to areas being pushed out of an essential central functional regime of gamma power. This further suggests the importance of careful titration when using schizophrenia medication; it would be easy to overshoot a target gamma-power zone.

Methods

Model

The full model is available on modelDB#258738 as http://modeldb.yale.edu/258738. It is based on refs. [43,78], modelDB#139421. The model was implemented in NEURON 7.4[79] running in parallel on eight 2.67 GHz Intel Xeon quad-core CPUs[80]. Result robustness was tested by using multiple random wirings and random stimulation input patterns for each parameter set, resulting in around 1900 simulations.

The model consisted of 800 pyramidal (PYR) neurons, 200 oriens-lacunosum moleculare (OLM) interneurons, and 200 parvalbumin-positive (PV) fast-spiking basket interneurons (Fig. 8). We used equal numbers of PV and OLM interneurons as PV interneurons provide the most prevalent feedback inhibition to pyramidal soma, while OLM interneurons provide the most dominant feedback inhibition targeting distal dendrites (reviewed in chapter 2 of ref. [109], and similar to other models[43,78,110]). These two interneurons are also hypothesized to generate theta-modulated gamma oscillations[110-112]. All neurons were modeled using multiple channels defined using Hodgkin–Huxley parameterizations. PYRs had five compartments comprising soma, apical dendrite, and basal dendrite. OLM and PV cells were single compartment. All neurons contained leak, fast sodium, delayed rectifier, and Ih. PYRs also had A-current, and had increasing Ih conductance up apical dendrite[113]. OLM added Ca2+-gated K2 current and high-threshold Ca2+ current with intracellular calcium dynamics.

Fig. 8

Schematic diagram for the in silico CA3 network model showing targets for multi-drug target manipulations.

Population of pyramidal neurons (n = 800) is represented by the red cell, PV basket interneurons (n = 200) by the green cell, and OLM interneurons (n = 200) by the blue cell. Dotted lines represent the random input driving the pyramidal neuron population (driving inputs). Output represents the spiking output of pyramidal neurons. MS is medial septum, providing inhibitory input onto the PV basket and OLM interneurons. Numbers on connections are convergence ratios for connectivity between different populations. Cyan boxes show the molecular targets that are being investigated in this study and their locations: on OLM interneurons, on pyramidal neurons, and on pyramidal neurons and PV basket interneurons. Activation and inactivation time constants for each of the receptor types are shown.

PYRs projected both AMPARs and NMDARs on all cell types, with mid-apical projections onto other PYRs. PVs and OLMs projected to GABAAR: PV → PYR soma; OLM → PYR distal apical. OLM, PV received GABAAR input from medial septum (MS) at 6.7 Hz (theta). Synapses were modeled as double exponential mechanisms with parameters from ref. [110] (Table 1). Background activity was simulated low amplitude Poisson input to all populations[114] (Table 2). Five seconds of simulation time were used for analysis, after discarding the first two seconds to allow network dynamics to stabilize. Local field potential (LFP) was calculated as voltage difference between apical and basal dendrites of all PYRs.

Table 1

Synaptic parameters for neuronal connectivity.

Presynaptic	Postsynaptic	Receptor	τ1 (ms)	τ2 (ms)	Conductance (nS)	Convergence
Pyramidal	Pyramidal	AMPA	0.05	5.3	0.02	25
Pyramidal	Pyramidal	NMDA	15	150	0.004	25
Pyramidal	Basket	AMPA	0.05	5.3	0.36	100
Pyramidal	Basket	NMDA	15	150	1.38	100
Pyramidal	OLM	AMPA	0.05	5.3	0.36	10
Pyramidal	OLM	NMDA	15	150	0.7	10
Basket	Pyramidal	GABAA	0.07	9.1	0.72	50
Basket	Basket	GABAA	0.07	9.1	4.5	60
Basket	OLM	GABAA	0.07	9.1	0.023	1
OLM	Pyramidal	GABAA	0.2	20	72	20

Table 2

Synaptic parameters for modeling background activity.

Cell	Section	Synapse	τ1 (ms)	τ2 (ms)	Conductance (nS)
Pyramidal	Soma	AMPA	0.05	5.3	0.05
Pyramidal	Soma	GABAA	0.07	9.1	0.012
Pyramidal	Dend	AMPA	0.05	5.3	0.05
Pyramidal	Dend	NMDA	15	150	6.5
Pyramidal	Dend	GABAA	0.07	9.1	0.012
Basket	Soma	AMPA	0.05	5.3	0.02
Basket	Soma	GABAA	0.07	9.1	0.2
OLM	Soma	AMPA	0.05	5.3	0.0625
OLM	Soma	GABAA	0.07	9.1	0.2

Manipulations

Target locations for NMDAR and Ih pathological manipulations were based on our prior work: NMDAR on OLM: modelDB#139421[34,78]; Ih on PYR, PV modelDB#151282[43]. We manipulated GABAARs on PYRs based on localization evidence from mRNA study in schizophrenia[52] (Table 3).

Table 3

Manipualtions done in the study.

Channel	Cell	Manipulation
NMDAR	OLM	⇑
GABAA	Pyramidal	⇕
Ih	Pyramidal & Basket	⇕

Output measures

LFP power spectrum density (PSD) was calculated using Welch method (Python Scipy signal module)[115] after removing DC component. Gamma power was measured as a 25–50 Hz band. Normalized transfer entropy (nTE) was used to measure information flow with 15 ms binning[88,116]; see ModelDB#136095 and Fig. 9. We measured phase-amplitude coupling between theta phase and gamma amplitude using modulation index (MI)[117,118]. Firing synchrony was measured using population coefficient of variance (popCV): for interspike intervals (isi) of the neuronal populations[81].

Fig. 9

Schematic diagram for steps to calculate transfer entropy (TE) from input to output spikes of a single pyramidal neuron.

a Input signal is represented by red arrows (darker red for Past, lighter red for Future). Output signal is represented by green arrows (darker green for Past, lighter green for Future). Spikes in both input and output signals were binned, and the spike counts were used to generate probability distributions to calculate entropy. b The entropies of the signals are overlapped, represented by areas of overlap of the sets in the Venn diagram. The overlap of the entropies is partly because of the flow of information between the different signals, and partly because of chance. c The blue region, representing the overlap between H(inputPast) and , is the mutual information (MI) between the Past of the input signal and the Future of the output signal (). The striped portion of the blue region represents TE from the Past of the input signal to the Future of the output signal. TE was calculated as the MI between H(inputPast) and , given the Past of the output signal ().