Alterations in Schizophrenia-Associated Genes Can Lead to Increased Power in Delta Oscillations.

Genome-wide association studies have implicated many ion channels in schizophrenia pathophysiology. Although the functions of these channels are relatively well characterized by single-cell studies, the contributions of common variation in these channels to neurophysiological biomarkers and symptoms of schizophrenia remain elusive. Here, using computational modeling, we show that a common biomarker of schizophrenia, namely, an increase in delta-oscillation power, may be a direct consequence of altered expression or kinetics of voltage-gated ion channels or calcium transporters. Our model of a circuit of layer V pyramidal cells highlights multiple types of schizophrenia-related variants that contribute to altered dynamics in the delta-frequency band. Moreover, our model predicts that the same membrane mechanisms that increase the layer V pyramidal cell network gain and response to delta-frequency oscillations may also cause a deficit in a single-cell correlate of the prepulse inhibition, which is a behavioral biomarker highly associated with schizophrenia.

Introduction

A recent genome-wide association study (GWAS) of schizophrenia (SCZ) identified 145 loci exceeding genome-wide significance (Pardiñas et al. 2018), several of which were missed by the previous consortium GWAS (Ripke et al. 2014). At least ten of these loci implicate genes that encode ion channels and calcium transporters and are thus major contributors to the electrical properties of neurons (Ripke et al. 2014; Pardiñas et al. 2018). The biological mechanisms by which the disorder symptoms emerge are unknown, but there is a growing body of data, including recent genetic evidence (Devor et al. 2017), that SCZ is a disorder of synaptic plasticity and connectivity (Stephan et al. 2006) or cortical excitability at large (O’Donnell 2008; Hasan et al. 2013). This is reflected in clinical biomarkers and endophenotypes of SCZ that can be observed and quantified using electroencephalography (EEG). In parallel, biophysically detailed modeling has reached a level where the genetic contributions to distinct ion-channel species of the neurons can be addressed and the neuron’s signature in macroscopic EEG signals can be predicted (Grillner 2014; Wang and Krystal 2014). In the present modeling work, we combine these two advances and show how common variants of ion-channel and calcium transporter-encoding genes may affect the macroscopic electrical signals, leading to a frequently observed clinical phenotype in SCZ.

The biomarkers and endophenotypes of SCZ include deficits in prepulse inhibition (Turetsky et al. 2007) and mismatch negativity (Umbricht et al. 2003) and modified power spectra of delta (Sponheim et al. 1994; Duan et al. 2015) (0.5–5 Hz), alpha-theta (Sponheim et al. 1994; Hong et al. 2008) (4–13 Hz) and gamma (Hall et al. 2011) (25–70 Hz) oscillations. The prepulse inhibition of the startle reflex is defined by a decrease of the amplitude of the startle reflex when a startling stimulus is preceded 30–300 ms by a weak prestimulus (Turetsky et al. 2007). SCZ patients and their first-degree relatives exhibit a robust deficit in this inhibition, and due to the similarity of the startle pathway across mammalian species, this phenotype has also been intensively studied using animal models (Turetsky et al. 2007). The increase in delta power, in turn, is unique (Michie et al. 2002; Braff et al. 2008; Hong et al. 2008; Hall et al. 2011; Duan et al. 2015) in that it has been found only in SCZ patients and not in their first-degree relatives and has thus been hypothesized to be causal (Duan et al. 2015; Lisman 2016). In general, delta oscillations are prominent during sleep and mental tasks. Supporting the prominent role of delta-oscillation alterations in SCZ pathology, patients with SCZ have repeatedly been shown both to suffer from sleep abnormalities (Keshavan et al. 1998) and perform worse than controls in many cognitive tasks, such as arithmetic tasks, working-memory tasks and tasks measuring cognitive control (reviewed in Hinkley et al. 2010 and Schwarz et al. 2016; cf. Duan et al. 2015).

The details of the mechanisms behind delta oscillations remain to be elucidated, but they seem to have two components generated in the thalamus and neocortex (Neske 2015). The thalamically generated delta oscillations stem largely from the intrinsic properties of the thalamocortical neurons. The cortically generated delta oscillations likely rely on the properties of layer V pyramidal cells (L5PCs) (Carracedo et al. 2013; Neske 2015; Blaeser et al. 2017). In vivo, L5PCs show transitions between “up” and “down” states, which correspond to states of high and low excitability (and firing rate), respectively (Harmony 2013). These transitions occur in slow (<1 Hz) and delta frequencies, giving rise to the amplification of these frequencies in EEG power spectra (Neske 2015). The down state has been proposed to be caused, at least in part, by L5PCs entering the afterhyperpolarization (AHP) period after synchronized bursting during the up state (Harmony 2013; Neske 2015). This hypothesis underlines the possibility that the properties of delta oscillations are largely dependent on not only synaptic connections but also non-synaptic ion channels contributing to the AHP period. If true, this would strengthen the link between altered delta oscillations in SCZ and the recent discoveries of gene variants related to ion-channel and Ca2+ transporters associated with SCZ (Ripke et al. 2014; Devor et al. 2017). 
In this work, we use established models of L5PC networks to show that the excitability of the networks and the magnitudes of their responses to inputs that oscillate in a delta frequency are significantly altered by modeled variants of SCZ-associated, ion-channel and Ca2+ transporter-encoding genes. We do this by incorporating the previously constructed (Mäki-Marttunen et al. 2016) models of SCZ-associated variants in a circuit model of thick-tufted L5PCs (Hay and Segev 2015), using reduced-morphology representations of the L5PCs to speed up the simulations (Mäki-Marttunen et al. 2018). Due to the lack of functional genomics data on common single-nucleotide polymorphisms (SNPs), we do not know which of our model variants represent cellular-level phenotypes of the SCZ-associated SNPs—however, our findings based on a large (N = 101) set of model variants suggest a correlation between increased delta-oscillation power and a single-cell correlate of the deficit in prepulse inhibition. We complement these analyses by computational experiments where we study, based on our in-house blood sample data, the effects of altered expression of specific ion channels or calcium transporters on network dynamics. Importantly, we are able to estimate the EEG signature of our L5PC population using the model of Næss et al. (2017) and show that the effects of the SCZ-associated variants may be directly reflected in the EEG signal. Our approach thus bridges the gap between the levels of individual genes and macroscopic electrophysiological signals, proposing a novel mechanistic link between the newly identified risk genes and the clinical phenotype of increased delta-oscillation power.

Methods

L5PC Network Model

We employed the single-cell “Hay” model of thick-tufted L5PCs (Hay and Segev 2015) as the basis of our study. This model, built using extensive electrophysiological data from rat neocortex, includes a reconstructed morphology and descriptions of eleven types of ion channels as well as a simple description of the intracellular Ca2+ dynamics (Hay and Segev 2015). The model thus represents a medium- to high-complexity neuron model that is well suited for studying contributions of different ion channels to neural responses. In this work, we coupled this model with human in vitro electrophysiological data on ion-channel behavior from functional genomics literature, following the approach of Mäki-Marttunen et al. (2016). We used the circuit model consisting of 150 L5PCs with identical morphology, as presented in Hay and Segev (2015), with the following modifications. To reduce overall simulation times, we replaced the original Hay model with 196 compartments by a four-compartment neuron model. This reduced-morphology model was presented in Mäki-Marttunen et al. (2018), where the ion-channel conductances, passive membrane parameters and parameters governing the Ca2+ dynamics were fitted in a stepwise manner to data obtained from simulations of the Hay model. In addition to these changes, we corrected an error in the original model (Hay and Segev 2015) that was causing depletion of the pre-synaptic vesicles even when no release occurred. The model L5PC received AMPA, NMDA and GABA receptor-mediated background synaptic inputs and AMPA and NMDA receptor-mediated L5PC-to-L5PC synaptic currents. As in all Hodgkin–Huxley-based systems, the integration of these inputs will be affected by the ion-channel mechanisms—the effects of the model variants we use stem from alterations of this ion-channel-contributed integration. The single-cell and network models are presented in detail in Supplementary Sections 1.1.1–1.1.4. The NEURON software (Hines and Carnevale 1997) was used for simulating the model.

To confirm that our results are not specific to networks consisting of only excitatory neurons, we explored the effects of our model variants in networks where the L5PC population is randomly connected to an inhibitory basket cell population (N = 50). For the basket cells, we used the single-compartment model for fast-spiking interneurons (Pospischil et al. 2008), which were connected to each other and the L5PCs with chemical synapses, obeying the connectivity statistics from Markram et al. (2015). Furthermore, we connected the basket cells to each other with gap junctions, as suggested by experimental data (Galarreta and Hestrin 2002) (see Fig. S1 for an illustration of the activity in these excitatory–inhibitory networks and Supplementary Section 1.1.5 for details on the construction of these networks). The models for the effects of genetic alterations are presented in Supplementary Section 1.2. The methods for sampling oscillatory Poisson processes (needed in simulation of the responses of the networks to oscillatory inputs) are described in Supplementary Section 1.3, and Supplementary Section 1.4 describes the methods for quantifying the oscillations in the spike train data. Finally, Supplementary Section 1.5 presents the scheme for forward modeling of the EEG signal.

Modeling Effects of Genetic Variants

We followed the procedure of Mäki-Marttunen et al. (2016) in modeling the effects of SNP-like genetic variants, see Supplementary Section 1.2 for details. In brief, we first chose the set of genes of interest based on the SNP-wise P-value data of Ripke et al. (2014) and searched the literature for variants of these genes and their effects on cell electrophysiology (Supplementary Section 1.2.1). We restricted our survey on ion-channel or Ca2+ transporter-encoding genes that are likely to be expressed in L5PCs (see below), and as a result, obtained models for variants (see Massa et al. 1995; Ji et al. 2000; Gomora et al. 2002; Lesso and Li 2003; Murbartián et al. 2004; Tang et al. 2004; Hohaus et al. 2005; Mantegazza et al. 2005; Ficarella et al. 2007; Ishii et al. 2007; Vanmolkot et al. 2007; Cestèle et al. 2008; Stary et al. 2008; Cordeiro et al. 2009; Kudrnac et al. 2009; Link et al. 2009; Empson et al. 2010; Hu et al. 2010; Bocksteins et al. 2011; Depil et al. 2011; Giacomello et al. 2011; Tan et al. 2011; Volkers et al. 2011; Zhang et al. 2011; Fakira et al. 2012; Lieb et al. 2012; Azizan et al. 2013; Cestèle et al. 2013; Pinggera et al. 2015; Wemhöner et al. 2015) of the following genes: CACNA1C, CACNA1D, CACNB2, CACNA1I, ATP2A2, ATP2B2, SCN1A, HCN1 and KCNB1 (see Table S2). However, the effects of most literature-based model variants on the L5PC excitability were large (see Mäki-Marttunen et al. 2016, Figure 1 and Fig. S1), which is not expected from common variants. Therefore, following the procedure of Mäki-Marttunen et al. (2016), we imposed conditions on how much the neuron responses (to predefined stimuli) were allowed to vary and scaled down the parameter changes related to each of the literature-based model variants until these conditions were met. In practice, each change of parameter that was related to a considered literature-based model variant as presented in Table S2 was multiplied (on linear or logarithmic scale, depending on the type of the underlying parameter) by a scaling parameter c = c, where the parameter c was determined as the largest factor c [0, 2] for which the scaling conditions 1–5 (Supplementary Section 1.2.2) concerning single-neuron firing behavior were met. This gave us a set of 101 small-effect model variants (Table S3), which we then further downscaled with fixed factors ε and employed as models for effects of common variants of the underlying genes on L5PC electrophysiological properties. Throughout the work, we use the term “variant” to refer to a gene variant in the genome of the human or experimental animal and the term “model variant” to a model of a gene variant where certain parameters of the L5PC model are changed (as indicated above) to describe the effects of a gene variant.

Modeling Effects of Altered Gene Expression—Indications from Blood Sample Data

Many of the SCZ-related SNPs affect the levels of expression of the targeted proteins without affecting their structure (Ripke et al. 2014). In contrast to the models of common SNP-like variants that (in the case of ion channels) affected the voltage dependence and kinetics of the target channels, altered expression of an ion-channel-encoding gene could be modeled as a change in total maximum conductance of the underlying channel type. To do this, we needed data on which ion channels are likely to be under- or over-expressed in neurons of SCZ patients in comparison to those of healthy controls. Measuring the gene expression in human neurons can generally be carried out only in very restricted cases, but indicative information can be obtained from blood sample data (Sullivan et al. 2006). Here, we used data from the thematically organized psychosis (TOP) study sample (Dieset et al. 2012, 2015), Oslo University Hospital, including 338 SCZ patients (including schizoaffective and schizophreniform patients) and 263 healthy controls. A total of 61% of the patients (206/338) and 53% of the controls (98/184) were male. The patients (31 ± 10 years) were slightly younger than the controls (33 ± 10 years) (t-test, P-value = 0.034 < 0.05). Blood samples were collected using Tempus Blood RNA Tubes (Life Technologies Corporation, Carlsbad, CA, USA). Total RNA was extracted with ABI PRISM 6100 Nucleic Acid PrepStation (Life Technologies Corporation) and TEMPUS 12-port RNA Isolation Kit according to the manufacturer’s protocol.

We analyzed samples using Illumina Human Expression Arrays (Illumina HT−12) (Illumina, San Diego, CA, USA). We used multidimensional scaling and hierarchical clustering for regular quality control and removed outliers and multiple batch effects. After this, we performed a log2 transformation. We ignored probes that showed zero expression in more than 90% of the samples. We refer to our earlier work (Akkouh et al. 2018) for more details on the RNA microarray analysis and quality control. We fitted a linear model to the expression level with respect to age and gender and tested the difference of the residuals in healthy control versus SCZ patient data (t-test, P-value 0.05). The use of antipsychotics was not controlled. Of the 267 assays for genes encoding for ion-channels and Ca2+-transporters, 30 assays had significantly different expression in SCZ patients versus healthy controls. From these 30 assays, we concentrated on probes targeting genes that are expressed in the cortex and whose effects can be studied using our L5PC model: these were ATP2A2, ATP2A3, ATP2B1, CACNA2D2, CACNA2D3, CACNA2D4 and CACNB3. However, genes ATP2A3, CACNA2D2 and CACNA2D4 do not seem to be expressed in L5PCs (see below); thus, we disregard them in this work.

The effects of altered expression of these genes can be directly or indirectly represented in our L5PC model. The β subunit encoded by CACNB3 and the α2δ subunit encoded by gene CACNA2D3 can modulate the number of high-voltage-activated (HVA) Ca2+ channels inserted to the membrane as well as the conductance and kinetics of the channel (Dolphin 2003, 2013). Here, we modeled the effects of these genes by altering the maximal conductance of the ICaHVA current (over-expression ↔ increase in maximal conductance; Canti et al. 2001; Dolphin 2013). The gene ATP2A2 encodes for a subunit of the SERCA protein, which pumps free intracellular Ca2+ into the endoplasmic reticulum (ER) or sarcoplasmic reticulum (SR). As the description of the Ca2+ stores was not included in the single-L5PC model, the effect of this gene was modeled indirectly. The model includes a parameter describing the fraction γ of the incoming Ca2+ currents that remain in the sub-membrane area—the underlying process involves various Ca2+ buffers and transporters, but the number and efficiency of the SERCA pumps are likely to be an important factor to this quantity (Ji et al. 2000). Thus, we modeled the contribution of expression level of ATP2A2 by altering the value of this parameter (over-expression ↔ increase in γ). Finally, the gene ATP2B1 encodes a PMCA pump, which expels the Ca2+ ions from the intracellular medium to the extracellular matrix. The effect of altered numbers of these proteins was modeled as a changed value of the parameter τdecay that represents the time constant of the Ca2+ decay towards the resting-state value (over-expression ↔ decrease in τdecay; see Supplementary Section 1.1.2).

Of the genes of Table S2, ATP2A2 was the only SCZ-associated gene that also showed differential expression in our blood sample data. This observation is in line with a post-mortem study showing increased expression of ATP2A2 in hippocampi and prefrontal cortex of SCZ patients (Earls et al. 2012). Similar results have been previously obtained for a SCZ-associated variant rs1006737 in CACNA1C, but the direction of the effect was brain region specific (Bigos et al. 2010; Eckart et al. 2016), and thus, it is not surprising that no difference was found in the expression level in blood samples. Likewise, there are rare SCZ-associated variants of CACNA1I gene that have shown altered expression in embryonic kidney cells (Andrade et al. 2016); however, we did not find significantly different expression of CACNA1I gene in our SCZ versus control data (see Table 1).

Table 1

Expression of ion-channel or Ca2+ transporter-encoding genes.

Gene	RNA expression in HPA dataset (TPM)	Protein expression level in HPA dataset	Mouse L5 visual cortex RNA expression		P-value
			FPKM	TPM
(A)
ATP2A2	119.7	Medium	77.61 ± 41.53	125.08 ± 72.46	0.034 (↓)
ATP2A3	1.9	High	0.007 ± 0.28	0.011 ± 0.45	0.00061 (↓)
ATP2B1	143.6	High	58.86 ± 43.73	94.48 ± 71.12	0.050 (↓)
CACNA2D2	5.9	Low	1.37 ± 2.63	2.14 ± 4.12	0.00033 (↓)
CACNA2D3	16.2	N/A	44.66 ± 51.17	72.21 ± 85.80	0.019 (↓)
CACNA2D4	0.3	N/A	0.11 ± 1.18	0.18 ± 1.92	0.026 (↑)
CACNB3	33.6	N/A	18.91 ± 18.98	30.84 ± 32.43	0.010 (↓)
(B)
CACNA1C	8.3	Medium	1.85 ± 3.31	3.00 ± 5.54	0.17
CACNA1D	4.1	Medium	2.60 ± 4.30	4.26 ± 7.58	0.28
CACNB2	26.3	N/A	7.34 ± 10.59	11.82 ± 17.03	0.13
CACNA1I	4.7	N/A	0.87 ± 2.08	1.42 ± 3.52	0.32
ATP2B2	64.7	High	33.22 ± 21.20	53.75 ± 36.95	0.11
SCN1A	10.3	N/A	22.47 ± 19.49	35.63 ± 31.06	0.88
HCN1	4.9	Medium	10.71 ± 12.28	16.86 ± 19.28	0.39
KCNB1	21.6	Medium	9.95 ± 9.96	16.32 ± 18.56	0.66

The first column shows the (human) gene name. The second and third columns show the corresponding RNA and protein expression values in the Human Protein Atlas, respectively. The fourth and fifth columns show the RNA expression of the corresponding murine genes in cells whose somata are located in layer V of the mouse visual cortex. These data are from the Mouse Brain Atlas: the values in the fourth column are in fragments per kilobase million (FPKM), while the values in the fifth column are in transcripts per million (TPM). The sixth column shows the P-value from the t-test of differences between the gene expressions in healthy controls versus SCZ patients in our blood sample data: if there were multiple assays for a single gene, the one yielding the minimum P-value is shown. If the difference in expression between patients and healthy controls was significant (P-value 0.05), the arrow shows whether the gene was over-expressed (↑) or under-expressed (↓) in patients in comparison to healthy controls. (A) Data for the genes that showed significantly different expression in the blood when SCZ samples were compared with healthy controls. (B) Data for the genes identified by the SCZ GWAS data (see Table S2). One gene, namely, ATP2A2, overlapped in these sets—this gene is included in part (A).

Suggestions for SCZ-associated genetic effects on ICaHVA current are also found in the literature. In Maschietto et al. (2015), (cultured) human induced pluripotent stem cells (iPSCs) were found to differentially express CACNA1C (contributor to ICaHVA current) between healthy control and SCZ patient-derived cells. As our L5PC model does not differentiate between channels composed of different subunits, altered expression of CACNA1C would lead to parameter modifications qualitatively similar to that of CACNA2D3 and, thus, we do not explicitly consider the altered expression of CACNA1C.

Indications of Gene Expression in L5PCs

An important prerequisite for the eligibility of our approach is that the genes mentioned above be expressed in L5PCs. To date, there is no single study or database showing gene expression levels of all above genes in L5PCs specifically, but expression data can be combined from multiple sources to obtain a sufficient level of confirmation. Table 1A shows the expression values from the Human Protein Atlas for the genes indicated by the above blood sample analysis, and Table 1B for the corresponding data for genes of Table S2. However, data from the Human Protein Atlas only show average expression values across various cell types and should thus be supported by more specific expression data. The Mouse Brain Atlas database (another database by the Allen Brain Institute) provides layer-specific expression data from murine visual cortex. Out of the genes of Table 1, ATP2A3, CACNA2D2, CACNA2D4, CACNA1I and CACNA1C had low expression (<1 FPKM or <4 TPM) in layer V of visual cortex in the Mouse Brain Atlas database, while the rest of the genes had moderate or high gene expression values (fourth and fifth columns of Table 1). However, previous cell-type specific studies confirm the expression of some of these genes in L5PCs: mRNAs of CACNA1C and CACNA1I were found to be expressed in L5PCs of postnatal rats at different stages of development (Christophe et al. 2005). Furthermore, an early study showed the expression of CACNA1C in neurons from all layers of rat dorsal cerebral cortex (Hell et al. 1993). In conclusion, data from humans and mice indicate that all of the genes considered in the results section (CACNA1C, CACNA1D, CACNB2, CACNA2D3, CACNA1I, ATP2A2, ATP2B1, ATP2B2, SCN1A, HCN1 and KCNB1) are expressed in L5PCs.

Results

We performed simulations of 150 interconnected L5PCs, driving the network with randomly activated background synaptic currents. We repeated each simulation using different model variants of Table S2 (we also varied the scaling of each variant, see below). For each model variant, the properties of ion-channel activation of Ca2+ dynamics were slightly perturbed, leading to altered activation of the ion channels, and thus modified spiking and network dynamics. The results below show that many aspects of L5PC single-cell and network activity were affected by the SCZ-associated model variants.

SCZ-Associated Variants Can Increase the L5PC Network Excitability

The neurons were activated by AMPA, NMDA and GABA receptor-mediated background synaptic currents, whose activation times obeyed stationary, independent Poisson statistics as in Hay and Segev (2015). We varied the activation rate of the background synapses to analyze the gain of the network. The firing rate of the network depended both on the activation rate of the background synapses (Fig. 1A) and on the presence and type of the model variant altering the ion-channel or Ca2+ transporter functions (Fig. 1B). Most of the model variants affected the gain of the network, that is, the input/output relationship with respect to the activation rate of the background synaptic inputs (Fig. 1C). Extensive simulations revealed that 39 of the ε = 1/2 model variants (101 in total) of Table S3 (purple) fired on average in a higher frequency and 62 in a smaller frequency than the control network (Fig. S3). In Figure 1C, a model variant that produced the highest average firing frequency (see Fig. S3) for the ε = 1/2 scaling was chosen for each gene (except for ATP2B2 for which all positively scaled variants had a weaker gain than the control). As expected, a combination of these model variants radically increased the network gain (the median of the spike frequency distribution at input rate r = 1.0 was significantly different from that in the control network, U-test, P < 0.01, 10 vs. 4 samples), despite the small effects of the individual model variants (Fig. 1C). The effects of the model variants on network gain were in line with their effects on single-cell excitability (insets of Fig. 1C).

Figure 1.

Model variants affect network excitability. (A) Population spike trains for control networks with different event rates of background synaptic inputs. The black spike train shows the response of the model to the inputs used in the original model (Hay and Segev 2015), that is, when excitatory and inhibitory synapses were activated on average at frequencies f = 0.72 and 7.0 Hz, respectively (the distribution of the inter-event intervals was stationary Poisson). In the other examples, these rates were multiplied by the shown factor r. (B) Population firing rate of networks with r = 1.0 as a function of time, smoothened using a 25-ms Gaussian kernel. Blue curve shows the control network data, while the dashed magenta curve shows data from the CACNA1C model variant with ε = 1/2 scaling (corresponds to the CACNA1C variant of panel C). (C) Gain curves of variant networks, where the mean neuron-wise firing rate during an 11-s simulation is plotted against the rate factor r. Averaged over 3–10 samples. C, inset: f–I curves of single-L5PCs with the corresponding model variants (Mäki-Marttunen et al. 2016). These data were obtained by stimulating the neuron with a somatic DC with the amplitude shown on the x-axis—the y-axis shows the resulting spiking frequency (in Hz). Blue: control network (no variants), other colors: downscaled model variants with different downscaling parameters ε. The parameter changes related to each model variant are stated in Table S3.

Response to Inputs Oscillating in Delta Frequency and Emergence of Delta Oscillations in L5PCs Are Altered by SCZ-Associated Model Variants

Here, we show that the model variants affect the way inputs oscillating within the delta-frequency band are integrated by the network. To do this, we made the Poisson process event rate λ in the generation of background synaptic spike trains vary in time in a sinusoidal form with a 25% amplitude (see Supplementary Section 1.3 for details and Fig. S2 for an illustration). We altered the background oscillation frequency from delta to beta range to analyze the resonance of the network to different frequencies. Different baseline frequencies led to different L5PC network firing responses (Fig. 2A). The response of the network was quantified using standard Fourier spectral analysis (Supplementary Section 1.4). Typically, the model variants affected the amplitudes of the power spectra of the population spike trains (Fig. 2B). The control network was previously shown to resonate strongly with inputs oscillating in a frequency of 0.7–2.5 Hz (Mäki-Marttunen et al. 2018). The L5PC networks expressing the SCZ-associated model variants (the same variants as in Fig. 1C) preserved this feature but, except from the network with KCNB1 variants that had negligible effects, altered the amplitude of the response (Fig. 2C). The median of the power spectrum amplitude at baseline frequency 1.5 Hz in the network with combination of model variants (ε = ) was significantly different from that in the control network (U-test, P < 0.01, 8 vs. 11 samples). The impacts of CACNA1I, SCN1A and HCN1 model variants were small but non-negligible (Figs 1C and 2C)—stronger versions of these variants (i.e., corresponding variants with larger scaling parameters) showed more pronounced effects on both L5PC network gain and response to delta oscillations (Fig. S4). The results for the combination of model variants were replicated using the full-morphology model of Hay et al. (2011)—Fig. S5 shows these results and confirms that the two models gave qualitatively similar results.

Figure 2.

Model variants affect the magnitude of the response of the network to oscillatory inputs. (A) Population spike trains for control networks where the event times of the background synaptic inputs were non-stationary (oscillatory) Poisson processes. The shaded areas show the times at which the λ of the oscillatory Poisson process was larger than average. The black spike train represents the L5PC network response to the input oscillation frequency 1.5 Hz, which caused the largest response in the network. (B) Fourier power spectra of the CACNA1C model variant (purple) and control (blue) network spike train with the input oscillation frequency 1.5 Hz. Only spikes from 2000 ms onward were considered when calculating these spectra. Inset: zoomed-in view on the peak at the frequency corresponding to the input oscillation. (C) The power of the frequency component corresponding to the input oscillation frequency (e.g., the amplitude of the peak in the inset of panel B) plotted against the input oscillation frequency. Averaged over 6–13 (mean 11.9) samples. Background oscillation frequencies f = 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 7.5, 10.0 and 15.0 Hz were considered. Blue: control network (no variants), other colors: downscaled model variants with different downscaling parameters ε.

In addition to responding strongly to delta-frequency inputs, the L5PC network is also hypothesized to have a key role in producing the cortical delta oscillations (Neske 2015). In Mäki-Marttunen et al. (2018), we demonstrated that given stationary (non-oscillatory) Poisson inputs, the L5PC network firing rate showed oscillations in the delta-frequency range. These oscillations were robust and depended on the SK currents (Mäki-Marttunen et al. 2018). In this work, we showed that the amplitude of the delta-frequency oscillations originating in the L5PC network was significantly altered by the combination of SCZ-associated model variants (Fig. S6A). In addition to the continuous-time continuous-frequency Fourier transform applied in Fig. S6A for determining the power spectra (see Supplementary Section 1.4), these oscillations could be observed using a range of FFT-based methods (Fig. S6B-E).

Effects of the SCZ-Associated Model Variants Are Observed in the Predicted EEG Signals

We predicted the EEG signal measured on the scalp (12.5 mm above the L5PC somata) using the head model consisting of four concentric spheres, where each layer had a constant electrical conductivity (Srinivasan et al. 1998). In this approach, the membrane currents of each compartment in each neuron were measured in order to calculate the dipole moment of the network and the EEG signal was evaluated based on these dipole moments using an analytically tractable model (Næss et al. 2017) (see Supplementary Section 1.5 for more details). In single-cell simulations, the events leading to somatic action potentials (Fig. 3A) and amplified synaptic currents (Fig. 3B) were reflected both on the dipole moment measurements (Fig. 3C) and the EEG signature of the neuron (Fig. 3D). In network simulations, by contrast, the large number of action potentials made these signals noisy, but periods characterized by high firing activity (Fig. 3E) and thus stronger synaptic currents (Fig. 3F) could still be recognized from the dipole moment and EEG signals (Fig. 3G,H). We calculated the spectral powers for different frequency components of the EEG signal and averaged across many (Nsamp = 140) repetitions, using both control networks and networks with the model-variant combinations of Figure 1C. The network data did not fit any single power-law behavior, but they showed approximate 1/f scalings with powers α = 1 for frequencies smaller than 20 Hz and α = 3 for frequencies greater than 20 Hz, in a fashion similar to that observed in experimental data from cat cortical EEG during waking (Bedard et al. 2006) (Fig. 3I). These values of α fall within the theoretically predicted range [1, 3] (Pettersen et al. 2014). We observed that the model variants that caused an increased firing frequency (purple, gray) in Figure 1C also caused an increase in the low-frequency powers—including delta-frequency range (Fig. 3J,K). This is an expected result, as increased spiking implies an increased magnitude of the dipole moments and the EEG signal magnitude is linearly dependent on the dipole magnitude. An increased magnitude of the EEG signal relates directly to an increased power at zero frequency (Supplementary eq. 5) and indirectly to powers at non-zero frequencies as well (Supplementary eq. 4). It is noteworthy that the shape of the power spectrum of the EEG signal (Fig. 3J) shows a less steep decay with frequency (steady decrease from 2 to ca. 40 Hz) than that of the population spike train (Fig. S6, from 2 to ca. 10 Hz), although both data are based on the same network simulations. This low-pass filter effect may be explained by intrinsic dendritic filtering (Pettersen and Einevoll 2008; Lindén et al., 2010; Leski et al. 2013).

Figure 3.

Predicted dipole moments and EEG signals from single-cell and network simulations and spectral EEG power in control and model-variant networks. (A) Somatic membrane potential time series of a single L5PC receiving the background synaptic inputs but no L5PC-to-L5PC inputs. The red rectangle is zoomed in on in the inset. (B) Synaptic currents summed across the dendritic sections in the L5PC of panel (A). During a burst of action potentials, the dendritic membrane potential is elevated, leading to magnified GABAergic (positive) background synaptic currents (inset). (C) Time series of the y-component (upwards, toward the scalp) of the dipole moment obtained from the L5PC of panel (A). (D) EEG time series (electrode placed at the top of the scalp) from the single-L5PC activity of panel (A). (E) Population spike train of 150 interconnected L5PCs. (F) Synaptic currents summed across the L5PC population of panel (E). (G) Time series of the y-component of the dipole moment obtained from the network of panel (E). (H) EEG time series from the L5PC network activity of panel (E). (I) Spectral power (see Supplementary Eqs. S2 and S4–S5) of the EEG signal such as that in panel (F), starting from 2000 ms, averaged over Nsamp = 140 networks with different random number seeds. The black lines represent 1/f and 1/f3 power scalings. (J) The data of panel (I) reproduced on a semi-logarithmic scale. The blue curve shows the power spectrum of the control network data of panel (I), and other colors represent the data from the same combinations of model variants as in Figs 1C and 2C. (K) Zoomed-in view on the data of panel (J).

Increased L5PC Network Gain and Delta Power May Be Caused by Altered Expression of Ca2+ Channels and Transporters

A parametric analysis showed that altering the maximal conductance of HVA channels or changing the ratio in which Ca2+ ions were included in the sub-membrane volume (γ) or the decay rate of intracellular Ca2+ concentration (τdecay) significantly affected the spiking behavior of the network (Fig. S7). In the analyses above, some of these parameters were affected by variants that modeled, based on data from functional genomics (see Table S2), the effects of altered amino acid sequences of the underlying proteins. Nevertheless, the values of these parameters could also be altered to more directly represent increased or decreased numbers of the proteins. In our sample, SCZ patients had altered expression of seven genes controlling the Ca2+ currents and intracellular Ca2+ homeostasis, four of which are likely to be expressed in L5PCs (Table 1A). Here, we considered several single-parameter variants that were related to these genes and tested their effects on network dynamics. We found that a modest (10%) decrease or increase in parameters describing the HVA Ca2+ channel conductance (gCaHVA), the inclusion parameter γ, or the decay parameter τdecay notably affected the L5PC network gain and response to oscillations (Fig. 4). The change of HVA Ca2+-conductance parameter simulates the altered expression of α2δ-3 (encoded by CACNA2D3) or β3 (encoded by CACNB3) subunit, whereas the changes in parameters γ and τdecay simulate the altered expression of SERCA (ATP2A2) and PMCA (ATP2B1) pumps, respectively. The bold curves of Figure 4 show the direction of effect as suggested by the expression data (under-expression of CACNA2D3, CACNB3, ATP2A2 and ATP2B1 in SCZ): interestingly, all these effects except that of the under-expression of ATP2B1 corresponded to increased network gain and increased responses to delta oscillations.

Figure 4.

Single-parameter variants affect the network gain and responses to oscillatory inputs. The panels from left to right show the effects of the following parameters: gCaHVA (maximal conductance of HVA Ca2+ currents), γ (proportion of Ca2+ ions that enter into the sub-membrane domain of the total number of incoming Ca2+ ions) and τdecay (the decay time constant of Ca2+ concentration in the sub-membrane domain). The upper panels show the gain curves (see Fig. 1C), and the lower panels show the responses to oscillatory inputs (see Fig. 2C). Oscillatory inputs with frequencies f = 0.5, 0.625, 0.75, 0.875, 1.0, 1.25, 1.5, 1.75, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 7.5, 10.0 and 15.0 Hz were considered. The densely dashed curves show the effects of a 10% increase in the parameter value, as the sparsely dashed lines show the effects of a 10% decrease in the parameter value. The curves corresponding to the trend (under-expression) observed in the blood sample data are thicker than the curves corresponding to the opposite trend. The solid lines show the data from the control network simulations. The data were averaged over Nsamp = 5–6 (upper panel) or 12 (lower panel) network simulations with different random number seeds.

Increase in Delta Power Is Not Canceled by Counteracting Synaptic Scaling

The results shown thus far illustrated how model variants of the risk genes increased the overall L5PC network excitability, which (especially in the case of combined variants) led to an amplified response to oscillations across all tested frequencies (see Fig. 2). To show that the delta-band responses of the L5PC network, in contrast to higher frequency ranges, are highly susceptible to the effects of the model-variant combinations of Figures 1C and 2C, we counterbalanced the variant-induced increase in network excitability by reducing the numbers of thalamocortical and corticocortical background synapses to NsynE = 10 000α and NsynI = 25 000α (α < 1). Using the bisection method for each model variant separately, we searched for the coefficient α that made the network implemented with the variant combination produce as much spiking as the control network in the baseline conditions (namely, when background synaptic inputs followed stationary Poisson statistics with λglut = 0.72 Hz and λGABA = 7.0 Hz). This search resulted in coefficients α = 0.872 (ε = 1/2) and α = 0.929 (ε = 1/4). We showed that the variant-induced increase in delta-band power persists although the power of higher frequency ranges are brought back to control levels by the synaptic scaling (Fig. S8). In vivo, synaptic scaling of this magnitude could be easily induced during animal development by mechanisms for homeostatic plasticity, given that the target neurons were intrinsically—or for external reasons—too active (Turrigiano and Nelson 2004). To support the hypothesis that this may happen for thalamocortical inputs to L5PCs in SCZ, there is anatomical evidence that the fiber pathways of inputs from midline and anterior thalamic nuclei to prefrontal cortex are reduced in volume in SCZ patients (Lang et al. 2006; Lambe et al. 2007) and computational modeling evidence from neural mass based models fitted to magnetic resonance imaging data that indicated reduced connection strengths of thalamocortical connections in high-risk subjects with psychotic symptoms (Dauvermann et al. 2013).

We repeated the simulations of Fig. S8 using models where the L5PC population was coupled to an inhibitory population of size N = 50 (see Supplementary Section 1.1.5). These simulations confirmed that the increase in delta power is robust against synaptic scaling also in excitatory–inhibitory networks, both in the presence of weak (Fig. S9) and strong (Fig. S9) gap junctions.

Predicted Connection between Reduced Prepulse Inhibition and Elevated L5PC Network Excitability

The above analysis is concentrated on neuronal function and information processing in layer V of the cortex. Nevertheless, the effects of the model variants on neuron responses may be generalizable to many types of neurons due to the wide expression of the underlying genes in the brain (Table 1). Here, we considered the effects of our model variants on the neuron’s response to two successive stimuli in a setting where the response to the second, stronger stimulus may be inhibited by the response to the first stimulus. This “cortical single-cell prepulse inhibition” is a robust phenomenon that may employ neural mechanisms similar to those of the prepulse inhibition of the startle reflex (Mäki-Marttunen et al. 2016). The startle network underlying the behavioral prepulse inhibition is well characterized and converges to the caudal pontine reticular formation (PnC), which directly activates the motoneurons (Lingenhohl and Friauf 1994; Leumann et al. 2001; Bosch and Schmid 2006). The giant PnC neurons show a long (up to 140 ms) AHP following their activation by inputs from neurons in the pedunculopontine nucleus (PPN) (Homma et al. 2002), which are activated by acoustic inputs (Rohleder et al. 2014). However, lesions of PPN do not completely block prepulse inhibition (Koch 1999), which favors the hypothesis that non-synaptic, voltage-gated ion channels (such as the ones studied here) residing in PnC neurons can also be important contributors to the prepulse inhibition.

We used the single-cell model of Hay et al. (2011) with reconstructed morphology. The neuron was equipped with a combination of the model variants of Figures 1C and 2C. In a similar fashion to Mäki-Marttunen et al. (2016), we randomly distributed 3000 excitatory synapses across the apical dendrite of the model L5PC. These synapses were synchronously activated at time t = 0 and later at time t = tI—however, unlike in Mäki-Marttunen et al. (2016), we here used a subthreshold amplitude for the first stimulus. The combination of model variants that increased the network firing frequency (see Fig. 1C) decreased the threshold conductance needed for the second input to induce a somatic spike in a L5PC (Fig. 5A). This can be interpreted as a deficit in the prepulse inhibition of the neuron. Accordingly, the combination of model variants that decreased the network firing frequency increased the threshold conductance for the second stimulus (Fig. 5A). Similar results were found for the single-parameter variants representing altered gene expression (Fig. S10). These behaviors were observed for inter-stimulus intervals ranging from approximately 40 to 120 ms, while for larger intervals, the trend was non-existent or even the opposite (by “opposite” we mean stronger prepulse inhibition than control, not prepulse facilitation; Wynn et al. 2004). We confirmed this behavior by calculating the correlation coefficients between the two phenotypes across all model variants of Table S3: the threshold conductance for a second stimulus, given 60 ms after the first stimulus, was negatively correlated with both the average network firing frequency (correlation coefficient −0.671, P-value 1.5 × 10−14) and the amplitude of the 1.5-Hz power of the population spike train given an input oscillating in a 1.5-Hz frequency (correlation coefficient −0.674, P-value 1.1 × 10−14) (Fig. S11A,B). By contrast, the threshold conductance for a second stimulus given 300 ms after the first stimulus was weakly correlated with the two network measures (correlation coefficients with the average network firing frequency and response to delta oscillations were 0.270 and 0.296, respectively, and the P-values were 0.0063 and 0.0027) (Fig. S11).

Figure 5.

The single-L5PC model predicts that the combination of variants increasing the network firing causes deficits in prepulse inhibition. (A) The solid curves show the threshold conductance for the second input arriving simultaneously to the 3000 synapses located at the apical dendrite. The dashed lines show the threshold conductance for the same input at rest; thus, the solid curves converge toward the dashed lines. For the positively scaled model variants (ε = 1/2 and 1/4), a deficit in prepulse inhibition, that is, lowered threshold conductance for the second input, was observed when the interval between the first and the second inputs is approximately 40–120 ms. The insets show the somatic membrane potential traces following the two stimuli, given different intervals and synaptic conductances of the second stimulus. If the two inputs were close in time (ISI = 40 ms), the neuron always fired at least one action potential shortly after the second stimulus. Likewise, if the interval between the two inputs was long (ISI = 100 ms), the neuron fired when the conductance of the second input was either g = 0.1 or 0.15 nS (ε = 1/2 variant fired also for g = 0.05 nS). By contrast, there was a lot of variation between the model variants in what happens for intermediate interval (ISI = 60 ms): ε = -1/2 variant did not fire for any of the three amplitudes of the second stimulus, ε = 1/2 fired for both g = 0.1 and 0.15 nS, and the behaviors of the other variants fell between these two extremes. (B) The SK current, recorded at the “hot zone” of Ca2+ channels (720 μm from the soma at the apical dendrite) in the single-L5PC model of Hay et al. (2011), both in the control case and when the combinations of model variants were used. Only the prestimulus (90% of the threshold conductance) was applied here.

The ion channel that has the largest contribution to the prepulse inhibition in the Hay model is the SK channel. As the SK current is Ca2+ dependent, large SK currents follow the activation of voltage-gated Ca2+ channels, even in the absence of a somatic action potential (Fig. 5B). While other K+ channels also contribute to hyperpolarizing the neuron, especially during the first 80 ms following the prestimulus, only the SK conductance in the Hay model has both the long decay and large amplitude needed to effectively inhibit the subsequent stimulus (data not shown). As a final step of our study, we replaced the subthreshold prestimulus by a suprathreshold stimulus and quantified the strength of the prepulse inhibition in terms of threshold conductance of the second stimulus for eliciting a second spike, as done in Mäki-Marttunen et al. (2016). Interestingly, the shape of the prepulse inhibition curve and the effects of the model variants remained qualitatively the same when this suprathreshold prestimulus was used (Fig. S12). These observations highlight the important role of dendritic SK currents both when the neuron is near to firing a somatic spike and when crossing the spiking threshold.

Finally, we also analyzed how the observed effects of SCZ-associated model variants on delta powers and single-cell prepulse inhibition were different from effects of other model parameter modifications that alter single-L5PC excitability but cannot be attributed to the effects of SCZ-associated genes. Fig. S13 shows that alterations in the membrane capacitance (Cm), length-to-diameter ratio (L/d; keeping the membrane area of each compartment constant) and axial resistance (Ra) of the neurite compartments and the recovery time from short-term depression (τrec) of the background synapses cause alterations in the predicted network gain and delta-resonance powers. While our approach does not allow considering the effects of τrec on prepulse inhibition, the parameters Cm, L/d and Ra also affected the predicted strengths of the prepulse inhibition. Network excitability and delta-resonance powers could be increased by increasing τrec (Fig. S13A,B) or decreasing Cm (Fig. S13), while changing L/d (Fig. S13) and Ra (Fig. S13) had mixed effects. The increase in network excitability by an increase in τrec is explained by the larger effect of this parameter on inhibitory background synaptic transmission than on the excitatory one, as the inhibitory synapses were activated at almost a ten-fold rate compared to the excitatory synapses (7.0 vs. 0.72 Hz at baseline rate factor r = 1). The increase in network excitability by a decrease in Cm is an expected result, as smaller Cm generally causes more rapid firing. Decreasing Cm also seemed to induce a modest positive shift in the peak resonance frequency, and vice versa (Fig. S13), which is consistent with previous in silico observations on CA1 pyramidal cells (Booth et al. 2016). The predicted prepulse inhibition was weakened by a decrease in Cm (Fig. S13), Ra (Fig. S13) or L/d (Fig. S13). Taken together, only a decrease in Cm produced results qualitatively similar to those obtained with the combination of model variants (Figs 1C,2C and 5). However, we found that compensating the decrease in Cm by decreased numbers of background synaptic inputs (Fig. S13) abolished the effects on delta-resonance powers (Fig. S13), unlike in the case of SCZ-associated model-variant combination (Fig. S8).

Discussion

Using computational modeling, we showed that altered expression or properties of the SCZ-associated, voltage-gated ion-channel or Ca2+ transporter-encoding genes lead to increased network excitability (Fig. 1) and delta power (Figs 2 and 3), which is a phenotype frequently observed in patients with SCZ. Many of these effects are dependent on the activation of the small-conductance, Ca2+-activated K+ channels (SK channels), which are expressed in L5PCs (Rudolph and Thanawala 2015). While these observations were expected based on the previous modeling results showing that SCZ-associated variants may alter single-L5PC excitability (Mäki-Marttunen et al. 2016, 2017), we here showed that the changes in delta power are robust against a compensatory change in synaptic strengths (Fig. S8) and the presence of an inhibitory population (Fig. S9). Furthermore, we showed that the same model variants that increased the delta power in the L5PC network also hindered the inhibition of the second apical stimulus in single-cell simulations of an L5PC (Fig. 5), mimicking the deficit in prepulse inhibition that is often observed in SCZ.

Despite the recent advances in iPSC reprogramming (Falk et al. 2016) and deeper characterization of SCZ endophenotypes (Braff et al. 2008), we are lacking cellular-level phenotypes that on one hand could be mapped down to genetic level and on the other hand would closely correspond to symptoms or phenotypes at the brain level. In this work, we proposed the increased intrinsic L5PC activity as a candidate for such a phenotype. Altered L5PC activity has previously been proposed as a potential mechanism for hallucinations (Larkum 2013), which are a positive psychotic symptom although not all patients of SCZ experience them. Furthermore, morphological alterations have been observed in L5PC of the prefrontal cortex in SCZ (Broadbelt et al. 2002; Black et al. 2004). On a circuit level, L5PCs contribute to the generation of delta oscillations (Neske 2015) and, thus, an increased intrinsic L5PC activity would be expected to increase the delta power in EEG recordings. Indeed, clinical studies frequently report increased powers of delta-frequency oscillations in patients of SCZ when compared with healthy controls (Clementz et al. 1994; Duan et al. 2015), although opposite results have also been found (reviewed in Hunt et al. 2017) (see Supplementary Section 3 for an extended discussion on this topic).

Previous computational studies on the topic have been performed using less detailed biophysical modeling, such as dynamic causal modeling (DCM). In one study, patients with psychosis and healthy controls performed an oddball paradigm task during EEG (Díez et al. 2017). The researchers used DCM to construct a neural mass model that predicted a decrease in frontal inhibitory connections in patients with psychosis, which led to local hyperexcitability of superficial pyramidal cells. Using similar methods in a mismatch negativity paradigm, an increase in the excitability of superficial pyramidal cells of inferior frontal gyrus was predicted in patients with psychosis in Ranlund et al. (2016). However, due to the spatial overlap between deep layer and superficial layer neuron populations, their contributions to the EEG signal can easily be mixed, and thus, we suspect that the predictions obtained for superficial populations in Ranlund et al. (2016) and Díez et al. (2017) may apply to L5PC populations as well. Network models more detailed than that of Ranlund et al. (2016) and Díez et al. (2017) (based on the microcircuit model of Bastos et al. 2012) could be used to more accurately characterize the source of neural alterations observed in patients with psychosis, which would allow a more detailed comparison of the results of biophysically detailed modeling of L5PCs and DCM-based results. On the other hand, while our model considers many possible factors (expression and intrinsic properties of the ion channels, connectivities within and between the L5PC and inhibitory neuron populations) to the excitability of the population, the DCM-based approaches only considered the connectivities as plausible factors (Bastos et al. 2012; Ranlund et al. 2016; Díez et al. 2017). Therefore, co-presence of increased intrinsic excitability (such as that obtained by the ion-channel-encoding gene variants in this study) and increased inhibition may be completely unrecognized by the DCM. This highlights the importance of our biophysically detailed modeling approach to complement the DCM-based approaches.

A centerpiece of our modeling framework is the downscaling of literature-derived model variants—which typically had large phenotypic consequences—into models of common SNP-like variants. This was done in order to prevent large effects in the physiology of the studied cell (see Mäki-Marttunen et al. 2016). Due to the great number of identified SCZ-associated gene variants and their frequency in the healthy population (Ripke et al. 2014), it is generally assumed that a single variant alone does not cause large consequences for brain activity either on a cellular or on a behavioral level. Note, however, that rare SCZ-associated variants with significant effects (e.g., Andrade et al. 2016) exist as well. To complement these analyses based on functional genomics data, we showed similar results (increased L5PC network excitability and response to delta-frequency oscillation and decreased cortical single-cell prepulse inhibition) for single-parameter model variants inspired by gene expression data. These gene expression data were based on blood sample analyses of SCZ patients and healthy controls and, thus, do not accurately reflect the differences in gene expression in brain tissue as such. However, gene expression data in the blood has been shown to be at least moderately correlated with expression data in the brain and other tissues (Sullivan et al. 2006). Furthermore, quantifying the gene expression from blood samples largely avoids the problem of RNA degradation that biases the results obtained from post-mortem studies (Koppelkamm et al. 2011). We consider the results obtained using these two approaches together as a strong signal of a modified L5PC firing behavior in SCZ, but experimental genetic animal model and/or iPSC model studies may be needed to confirm this.

Our blood sample analysis suggested an under-expression of genes CACNA2D3, CACNB3, ATP2A2 and ATP2B1 in SCZ patients compared to healthy controls. While our model predicted that an under-expression of ATP2B1 weakens L5PC excitability, the model predicted that the under-expression of genes CACNA2D3, CACNB3 and ATP2A2 leads to increased L5PC gain and larger amplitude of delta-frequency responses. Our results thus suggest that altered expression of genes encoding subunits of Ca2+ channels and transporters could explain the observation of increased delta-oscillation amplitudes in SCZ. By contrast, our results based on the SNP-like model variants that were constructed using functional genomics and GWAS data do not tell anything about the direction of the effects—this is because we included both gain-of-function and loss-of-function variants and always considered both positively and negatively scaled model variants. Indeed, we did not optimize our variant models (other than the magnitude of the parameter change) to reproduce a certain cellular or network effect but aimed at an unbiased analysis by taking a large number of variants as given in the literature and using previously validated off-the-shelf neuron models. This choice reflects the current lack of knowledge on the GWAS-identified variants: apart from a few rare SCZ-associated variants, we do not know anything about their functional effects. However, our modeling approach could help in characterization of the effects of variants in different genes (see Fig. S4) and in finding correlations between different phenotypes (see Fig. S11). These types of multimodal findings could help in revealing the polygenic structure of SCZ by making genetically and mechanistically based links between phenotypes that were previously conceived as separate. What we mean by this is that biophysically detailed modeling could produce detailed information such as “variants in genes A, B and C lead to phenotypes P and Q through effects in cellular properties X and Y,” while purely statistical genetics approaches can only find genetic links (“variants in genes A, B and C increase the prevalence of phenotypes P and Q”) and modeling with less biophysical detail can only involve the more superficial levels of abstraction (“cellular properties X and Y cause phenotypes P and Q”). Nevertheless, being more data-oriented than biophysically detailed modeling, statistical genetics and higher-level modeling (such as DCM-based approaches) have an important role in validating or invalidating the predictions made by biophysically detailed modeling.

Limitations and Future Directions

Our predictions are based on the Hay model that accurately describes the electrophysiological features of the thick-tufted L5PCs (Hay et al. 2011). Although the parameters of the Hay model were not fitted to reproduce the neuron activity under different ion-channel blockers (Mäki-Marttunen et al. 2018), the ion-channel structure of the model is very similar to that in other models of L5PCs with only a few differences. The Almog model (Almog and Korngreen 2014) describes the L5PC behavior using a partly overlapping set of ion channels. The persistent K+ channels have a larger role in the Almog model, and although the model separately describes BK currents (another Ca2+-activated current species), the contribution of the SK currents to the neuron behavior is significant in the Almog model as well (Almog and Korngreen 2014; Mäki-Marttunen et al. 2017). The model of Papoutsi et al. (2013) does not specify the molecular background of its Ca2+-activated slow AHP current, but its predictions are in line with the predictions of the Hay model in that the blockade of voltage-gated Ca2+ channels increased the L5PC activity due to the consequent decrease in the Ca2+-activated AHP current (Papoutsi et al. 2013). Thus, although the molecular structure of the Ca2+-activated AHP current in L5PCs remains to be characterized in more detail, the modeling (and partly the experimental, see Papoutsi et al. 2013) literature is to a large extent in agreement in that these currents are significant in L5PCs, and that due to these currents, pharmacological or genetic alterations that reduce the Ca2+ currents (including blockade of Ca2+ channels) typically increase the firing activity of the L5PCs despite the fact that the Ca2+ currents themselves are clearly depolarizing rather than hyperpolarizing (Mäki-Marttunen et al. 2017).

In this work, we constrained our analysis on genes encoding subunits of voltage-gated ion channels and Ca2+ transporters, but future work should address the contributions of the SCZ-associated genes encoding receptors of neurotransmitters as well as plasticity-related intracellular signaling molecules (Devor et al. 2017). Our framework could also be integrated with computational studies of NMDA receptor hypofunction and GABA deficiency, which are among the mechanisms more traditionally considered in computational studies of pathophysiology of SCZ (Jadi et al. 2016; Krystal et al. 2017). A major challenge for future computational work is to include the immune pathways in the models of SCZ pathogenesis, as indicated by GWAS and post-mortem studies (Ripke et al. 2014; Van Kesteren et al. 2017).

In the prediction of the EEG signal (Fig. 3), the use of the simplified, linear geometry for the L5PCs may have resulted in over- or underestimated dipole moments of the neurons. Future studies should employ more detailed modeling of the EEG signals (Tveito et al. 2017; Hagen et al. 2018) and involve a more diverse population of cortical neurons (Markram et al. 2015) to bring forth more realistic network dynamics. More work is needed also to characterize how and under which conditions (resting-state or task-related EEG) different neural populations contribute to the EEG features typically quantified by experimentalists. However, our results form a proof of principle on how activity in L5PCs maps to the EEG signal and how the effects of genetic variants may be reflected in the signal.

In conclusion, our work is a proof of concept on how data from GWASs and functional genomics can be integrated with biophysical modeling to tackle challenging questions regarding the pathophysiology of SCZ. Novel experimental methods could be used to test and further refine our model predictions in vitro and in vivo. In particular, new genome editing tools and automated cell-patching methods (Kodandaramaiah et al. 2012) could allow efficient in vitro analysis of the electrophysiological properties of a large number of common SCZ-associated variants, and information from such an analysis could be used to develop new animal or iPSC models of SCZ. However, as SCZ is a massively polygenic disorder, we believe that innovative experimental approaches will be needed to pinpoint specific genetic mechanisms in contrast to other genetic mechanisms that have overlapping effects. Our biophysically detailed modeling approach shows great promise as a tool for uncovering polygenic cellular-level mechanisms, as it allows studying alterations of both one genetic pathway at a time and in combination with other pathways. Given enough biological detail in the model, compensatory mechanisms that take place in parallel with the studied genetic alterations can be implemented for a better description of the observed experimental phenomena, but unlike in most experimental approaches for polygenic analysis, these mechanisms can be fully controlled by the experimenter. Challenges for the field are more detailed representations of cellular and sub-cellular entities in the model as well as the adaptability of the models to larger spatial and temporal scales.

Click here for additional data file.