We employ two neuron models, the “Hay model” (Hay et al., 2011) and the “Almog model” (Almog and Korngreen, 2014). Both models are multi-compartmental Hodgkin-Huxley type of models with reconstructed morphologies from layer V thick-tufted pyramidal neurons. The Hay model includes the following trans-membrane currents: Fast inactivating Na⁺ current (I_Nat), persistent Na⁺ current (I_Nap), non-specific cation current (I_h), muscarinic K⁺ current (I_m), slow inactivating K⁺ current (I_Kp), fast inactivating K⁺ current (I_Kt), fast non-inactivating K⁺ current (I_Kv3.1), high-voltage-activated (HVA) Ca²⁺ current (I_CaHVA), low-voltage-activated (LVA) Ca²⁺ current (I_CaLVA), small-conductance Ca²⁺-activated K⁺ current (I_SK), and finally, the passive leak current (I_leak). The Almog model includes a slightly different set of trans-membrane currents: Fast inactivating Na⁺ current (I_Nat), non-specific cation current (I_h), slow inactivating K⁺ current (I_Kp), fast inactivating K⁺ current (I_Kt), HVA Ca²⁺ current (I_CaHVA), medium-voltage-activated Ca²⁺ current (I_CaMVA), small-conductance Ca²⁺-activated K⁺ current (I_SK), large-conductance voltage and Ca²⁺-gated K⁺ current (I_BK), and finally, the passive leak current (I_leak).

It is not perfectly clear which particular ion channel subunits underlie each of these currents. mRNAs of ion channel-encoding genes KCNA2, KCND2, KCND3, CACNA1A, CACNA1B, CACNA1C, CACNA1D, CACNA1E, CACNA1G, CACNA1H, CACNA1I, HCN1, and HCN2 were observed in L5PCs in a study of postnatal rat neocortices at various stages of development (Christophe et al., 2005). Of these, CACNA1A, CACNA1B, CACNA1C, and CACNA1D encode the alpha subunit of HVA Ca²⁺ channels and thus contribute to I_CaHVA, and CACNA1G, CACNA1H, CACNA1I contribute to I_CaLVA and possibly I_CaMVA¹, while the genes KCNA2, KCND2, and KCND3 encode K⁺ channel subunits that might contribute to the slow I_Kp current, and the genes HCN1 and HCN2 encode subunits for a non-specific ion channel (contributing to I_h). Expression of KCNC1 was observed in certain subpopulations of L5PCs (Akemann et al., 2004), and expression of KCNB1 and KCNB2 was observed in layer V or VI pyramidal neurons (Guan et al., 2007). Channels encoded by KCNB1 and KCNB2 subunits likely contribute to the I_Kp current due to their slow activation kinetics, while KCNC1-based channels form the I_Kv3.1 current included in the Hay model. Expression of muscarinic potassium channel-encoding genes KCNQ2 and KCNQ3 in L5PCs was reported in Battefeld et al. (2014), and these genes are known to contribute to the I_m current.

Expression of Na⁺ channel subunit-encoding genes SCN2A and SCN8A was observed in pyramidal cells of all layers of human epileptic tissue (Tian et al., 2014). In another study, expression of genes SCN1A, SCN2A, SCN3A, and SCN6A in L5PCs was observed (Whitaker et al., 2000), while expression of SCN1A was not observed in Tian et al. (2014). Of the genes encoding Na⁺ channel α subunits, SCN1A, SCN2A, and SCN3A, alongside with SCN9A, are tetrodotoxin-sensitive (Catterall et al., 2005) and thus form both the transient (I_Nat) and persistent (I_Nap) Na⁺ currents. Whether these genes contribute to the current I_Nat or I_Nap may depend on the modulatory subunits that they are associated with (Ma et al., 1997).

Expression of KCNN1 and KCNN2 has been observed in L5PCs, while expression of the third gene affecting the I_SK current, namely KCNN3, is weak throughout the neocortex (although more prominent in pyramidal neurons in hippocampal region CA1) (Ballesteros-Merino et al., 2014). The α subunits of BK channels, encoded by KCNMA1 and giving rise to the I_BK current in Almog model, were found to be expressed in L5PCs in Benhassine and Berger (2005).

Most of the above genes encode the α subunit of the underlying ion channel. Many ion channels incorporate also modulatory subunits, the presence of which may change the kinetics or voltage-dependence of the ion channel. As an example, the β2 subunit of Ca²⁺ channels, encoded by CACNB2, associates with L-type Ca²⁺ channels (where the α pore is encoded by CACNA1S, CACNA1C, CACNA1D, or CACNA1F).

In addition to describing the dynamics of these trans-membrane currents, the Hay and Almog models also describe the dynamics of the intracellular Ca²⁺ concentration, [Ca²⁺]_i. According to the models, [Ca²⁺]_i is increased by the currents conducted by Ca²⁺ channels and otherwise decays toward a resting-state level of [Ca²⁺]_i. The extrusion of Ca²⁺ is contributed by many intracellular molecules, but two types of proteins possess a key role, namely, sarco/endoplasmic reticulum Ca²⁺ ATPase (SERCA) and plasma membrane Ca²⁺ ATPase (PMCA). The SERCA proteins pump cytosolic Ca²⁺ into sarcoplasmic or endoplasmic reticulum, which later releases the excessive Ca²⁺ into the cytosol to strengthen the Ca²⁺ transients caused by opening of the voltage-gated Ca²⁺ channels. SERCA pumps are encoded by three genes, namely, ATP2A1, ATP2A2, and ATP2A3, of which ATP2A2 is widely expressed in the brain (Baba-Aissa et al., 1998). PMCA proteins pump intracellular Ca²⁺ into the extracellular medium, and are encoded by genes ATP2B1, ATP2B2, ATP2B3, and ATP2B4, all of which are widely expressed in the brain (Stahl et al., 1992, 1994).

The distribution of the listed channels across the dendritic tree vary in these two models. In the Hay model, the basal dendrites only contain the leak and I_h current, and the initial axon segments only contain the leak current. The soma, by contrast, expresses all current types except for I_m, and the apical dendrite expresses all currents except for I_Nap, I_Kp, and I_Kt. The channel conductances are constants across the type of segment (somatic, apical, basal, axonal) except for the Ca²⁺ currents and the I_h current. The I_CaHVA and I_CaLVA currents are piece-wise constant such that the conductances is 10-fold or 100-fold, respectively, in the “hot area” of Ca²⁺ channels, located 685–885 μm from the soma in the apical dendrite. The I_h conductance increases exponentially by distance in the apical dendrite toward the distal end, but it is constant in soma and basal dendrites. In the Almog model, all active currents are location-dependent, such that the maximal conductances of I_Nat, I_CaHVA, I_CaMVA, I_SK, and I_BK decrease linearly from soma toward a predefined distance along the apical dendrite and then level off to constant. By contrast, the maximal conductances of I_Kp and I_Kt decrease exponentially from the soma toward the end of the apical dendrite, while that of the I_h current increases as a sigmoid function toward the end of the apical dendrite. All maximal conductances are non-zero constants along the basal dendrites.

Importantly, the two models are also different in terms of the firing behavior. One of the most obvious differences is that the Hay-model neuron fires tonically in response to somatic DC, while the Almog model expresses a rhythmically bursting (“chattering”) firing pattern, where each burst consists of 3–5 spikes. Both types of behavior are typical in the thick-tufted L5PC population, although the chattering firing pattern is less frequently observed (Chagnac-Amitai et al., 1990; Mason and Larkman, 1990; Yang et al., 1996). Both models are capable of describing back-propagating action potential-induced Ca²⁺ spike firing (BAC firing), where a coincidence or near-coincidence of a back-propagating action potential and apical stimulation generates a dendritic Ca²⁺ spike that propagates to soma to induce further spiking. In L5PCs, BAC firing is an important property that is hypothesized to underlie the special role of the cells in coupling feed-forward and feedback information (Larkum, 2013).

The morphologies accompanying the two models also differ from one another. The Almog model morphology is larger than that of the Hay model, both in terms of total dendritic length (13.5 vs. 12.7 mm) and total area (0.059 vs. 0.031 mm²), but is described by a smaller number of neurite sections (153 vs. 196).

We perform some of our simulations in the presence of dense synaptic background firing as expected in vivo. We do this following the methods of Hay and Segev (2015): we distribute 10,000 glutamatergic (AMPA and NMDA) and 2,500 GABAergic synapses across the apical and basal dendritic trees. The maximal synaptic conductances are 0.4 nS for glutamatergic and 1.0 nS for GABAergic synapses, and the reversal potentials are 0 and –80 mV, respectively. The AMPA and GABA conductances have dual exponential shapes with rise times 0.3 ms (AMPA) and 1 ms (GABA) and decay times 3 ms (AMPA) and 20 ms (GABA). The NMDA conductance has likewise a dual exponential shape with rise time 2 ms and decay time 65 ms, but it is filtered by a sigmoid function that describes the dependence of the conductance on Mg²⁺ concentration and membrane potential to take into account the Mg²⁺ block in the NMDA receptors (Jahr and Stevens, 1990). The activation times of the synapses follow a Poisson process with rate 0.72 Hz (glutamatergic) or 7.0 Hz (GABAergic). The synapses are short-term depressing with recovery time D = 800 ms and single-release consumption rate U_se = 0.6 (glutamatergic) or U_se = 0.25 (GABAergic).

We study the same set of genes as in our previous work (Mäki-Marttunen et al., 2016). Several of these genes were found to contain single-nucleotide polymorphisms (SNPs) associated with a high risk of SCZ (p-value smaller than 3·10⁻⁸ in the data of Ripke et al., 2014), namely, CACNA1C, CACNB2, CACNA1I, ATP2A2, and HCN1. Using a more relaxed threshold (p-value smaller than 3·10⁻⁵), we extended this set by the genes CACNA1D, SCN1A, KCNB1, KCNMA1 and ATP2B2. Compared to our previous study (Mäki-Marttunen et al., 2016), we left out the gene SCN9A and KCNN3 due their low expression in L5PCs, but included KCNMA1, which encodes the α subunit for the BK channel. The corresponding current, I_BK was described in the Almog model but not in the Hay model. We also omitted the gene KCNS3 that was included in Mäki-Marttunen et al. (2016) due to its marginal effects: it encodes a modulatory subunit of the Kv2 channels, which themselves showed small contribution to the firing behavior according to the model of Mäki-Marttunen et al. (2016).

It should be noted that we used the SNPs reported in Ripke et al. (2014) only to identify the above SCZ-related genes, and due to lack of functional genomics data, we could not include the actual physiological effects of the SCZ-related SNPs in our simulation study. Instead, we searched the literature for functional genomic studies reporting the effects of any genetic variants of these genes—in the case of voltage-gated ion channel-encoding genes, these effects typically included altered voltage-dependence and kinetics. Table A1 in Supplementary Material lists all such studies that we found in the literature where the effects of a variant of one of the above genes was measured in a way that could be directly implemented as a parameter change in one or more of our models. Not all of these genes have been shown to be expressed in thick-tufted L5PCs in specific (see above), however, they all show expression in the cortex (Mäki-Marttunen et al., 2018a).

As SCZ is a polygenic disorder, it is likely that the risk of obtaining the disorder is not caused by any of the SCZ-related common sequence variants alone, due to the small effect size for each. However, when sufficiently many of them are present, the combined effect could alter the electrophysiological properties. Although rare variants exist and contribute to the risk of SCZ by ~1–2% (Singh et al., 2017), there are not sequence variants that would be either necessary or sufficient for the emergence of the disease. To take this into account, we used the downscaling approach of Mäki-Marttunen et al. (2016) to model the effects of common SCZ-associated variants in L5PCs. In short, we simulated the neuron model with the parameter changes representing the variants as reported in Table A1 (Supplementary Material). If the variant dramatically altered the neural response, the parameters were brought closer to the control values so that there were no large difference between the firing behaviors of the control neuron and the mutant with the downscaled variant. The downscaling criteria were set so that no variant should alone dramatically alter the threshold amplitudes for chosen stimuli to initiate an action potential nor radically change the steady-state firing behavior of the neuron.

The variants of Table A1 (Supplementary Material) were scaled down to make sure the following four conditions were met:

In case one or more of the conditions I–IV were violated, the change of the model parameters related to the variant were scaled down, all parameters in proportion, to a threshold fraction c < 1 of the original effect (the smallest fraction with which the violation was observed). There are no data on how large parameter changes can be expected from common variants, and indeed, the chosen threshold values (±10% and ±15%) were hand-picked as representative values for mild deviations from control behavior. For this reason, we consider variants with different scalings where the threshold effect factor c is multiplied with another factor ϵ < 1. In this work, we consider the scaling factor values ϵ=12 and 14, and we also display the effects of the corresponding opposite variants ϵ=-12 and -14. For the variants that did not violate the conditions I–IV, we sought for the threshold effect up to twice the original effect (c < 2), and if the variant obeyed conditions I–IV for all values c < 2, we considered the original variant the ϵ=12 variant and applied other scalings with respect to that. The scaling coefficients for the Hay and Almog models are given in Table A2 (Supplementary Material), and Table S2 complements these data with the scaling coefficients for the altered Hay and Almog models. For more details on the variant scaling, see our earlier work (Mäki-Marttunen et al., 2016; Mäki-Marttunen et al., 2017).

Table S1 shows the amplitudes A₁, A₂, A_3a, and A_3b for each of the applied models. In the altered Hay models, the BAC firing properties according to which (among others) the unaltered Hay model was fitted in Hay et al. (2011) were radically changed. As a somatic square pulse current could induce a single spike in the unaltered Hay model, suprathreshold somatic currents always induced a burst of two or more spikes in the altered Hay models. These bursts were mediated by the more prominent activation of dendritic Na⁺ and Ca²⁺ channels than in the unaltered Hay model, which is shown by the fact that the dendritic membrane potentials are highly elevated immediately after the first spike unlike in the unaltered Hay model, as shown in Table S1. Such behavior has been reported previously in certain subpopulations of L5PCs (Amitai et al., 1993), although not in the data to which the Hay model and Almog model were fitted. Therefore, in the analysis below we are careful about the results obtained using altered Hay models, especially Hay-A₂–A₆, while we consider results from the unaltered Hay model as well as the unaltered and altered Almog models more established.

We used the NEURON software with Python interface for simulating the model neurons (Hines and Carnevale, 1997). Our simulation scripts are publicly available in ModelDB, entry 249463 (http://modeldb.yale.edu/249463).

We investigated the steady-state behavior of the different model neurons when a direct current (DC) was applied to the soma. Figures 2A–D shows the time courses and f-I curves (firing frequency as a function of the DC amplitude) for one variant as predicted by different models. Figures 2E,F shows the summary statistics of these data, and Figure 2G shows corresponding Hay-model data for a wider set of representative variants. These variants show notable differences from the corresponding control neuron data [Figure 2E (top panel) and Figure 2G]. However, the differences between the variant and control neuron f-I curves are limited by the explicit constraint posed by the scaling condition IV (see section 2). Although this constraint was only applied to the positively scaled variants (ϵ=12 and ϵ=14), their opposite variants (ϵ=-12 and ϵ=-14) mostly show well-constrained f-I curves as well (i.e., their f-I curves are not much further from the control data than those of the positively scaled variants). This indicates that the systems behave in a relatively linear manner with respect to model parameter alterations of this magnitude: for a more rigorous analysis on this (see Figure S1).

Figures 2E,F shows that for the effects of the chosen CACNA1C variant were qualitatively same across the models. Figure S2 shows that the same applies to the rest of the six variants of Figure 2G. We confirm this result in Figure S3, where the averages of the f-I curves for all variants of Table A2 (Supplementary Material) are plotted, and in Table S3, which shows the correlations of these quantities (for ϵ=12 variants) between the models. Figure S3C illustrates how these correlations are determined. The differences between the predictions of different models are caused by several factors. Firstly, the relationships between the half-activation voltages and slopes of the activation and inactivation curves were different between the Hay and Almog models. Therefore, a variant that simultaneously changes many of these parameters may have qualitatively different effects in the Hay models compared to those in the Almog models, see a discussion on this topic in Mäki-Marttunen et al. (2017). Secondly, in a bursting neuron model, the f-I curve shows steep increases when the number of spikes per burst increases by one, while the non-bursting spiking models have smoother f-I curves– thus, the firing rate of a lower-gain variant may for certain values of the input-current amplitude exceed that of the control neuron (see, e.g., the negatively scaled SCN1A variants for the Almog-model neuron in Figure S2B).

We confirmed our results also using the Hay model with spontaneously activated excitatory and inhibitory synapses (as in Hay and Segev, 2015). Figure S4A shows the membrane potential traces of control and a CACNA1C variant and Figure S4B shows f-I curves corresponding to the variants of Figure 2 when this in vivo-like noisy input was applied to the Hay-model neuron. Figure S4C shows the collection of data from all variants (corresponding to the data of Figure S3) in the noisy Hay-model neuron. The variant effects on noisy f-I responses (Figures S4B,C) were qualitatively similar to those in the absence of noise (Figures S2A, S3A). We also analyzed the effects of decreased NMDA and GABA receptor conductances on the Hay-model neuron firing under spontaneously activated synapses (Figure S4D). Importantly, the effects of the model variants of voltage-gated ion-channel and Ca²⁺ transporter-encoding genes were different from the impacts of altered synaptic inputs: while many model variants affected the L5PC gain (Figure S4B), alterations in synaptic conductances made the neuron fire with approximately the same increment or decline in firing rate across the tested perisomatic current amplitudes (Figure S4D). Decreasing the NMDA and GABA receptor conductances simultaneously did not appear to change the firing behavior of the L5PC, but kept the steady-state firing close to that observed in the control case (Figure S4D, right panels). Taken together, these results suggest that the steady-state firing of L5PCs is robustly (in a consistent manner across many neuron models, and also in presence of noise) affected by variants of SCZ-associated voltage-gated ion channel and Ca²⁺-transporter encoding genes, and that these effects are dissimilar to those caused by synaptic scaling.

Cortical neurons in vivo express a bimodal activity pattern that oscillates between a regime of high (“up” state) and low (“down” state) spiking activity in a slow (<1 Hz) frequency (Steriade et al., 1993). These states are prominent during anesthesia and sleep (Destexhe et al., 2007), but occur also in unanesthetized animals during quiet wakefulness (Petersen et al., 2003; Luczak et al., 2007). Most studies on up and down states apply electrophysiological measurements, but they have also been observed in Ca²⁺ imaging studies (see, e.g., Kerr et al., 2005; Grienberger et al., 2012). The origins and detailed characteristics of the up and down states are disputed, but it is known that L5PCs play an important role in the generation or maintenance of these states (Neske, 2016). The increased excitability in the up vs. down state could be contributed by at least three mechanisms: (1) intrinsic excitability (modulation of ion channels), (2) coordination of synaptic inputs, and (3) ephaptic coupling (Wilson, 2008; Anastassiou et al., 2011). In Hay et al. (2011), a simplistic model of the up state, where the proximal apical dendrite was depolarized with a long sub-threshold square-pulse stimulus, was presented, while the resting state of the neuron was considered the down state. This simplistic model of the up state could represent the effects of any of the above three mechanisms. In our previous study (Mäki-Marttunen et al., 2016), we showed by following this simplistic description of up and down states that SNP-like model variants of SCZ-associated genes affected the properties of integration of somatic and apical stimuli in both regimes. Here, we extended these results by considering a wider set of L5PC models and stimulation protocols, as well as a more in-vivo-compatible representation of the up and down states.

Following Hay et al. (2011) and Mäki-Marttunen et al. (2016), we implemented the simplistic up-state model such that the apical dendrite (at a distance of 200 μm from the soma) was stimulated with a 600-ms square pulse current (200 ms in Hay et al., 2011; Mäki-Marttunen et al., 2016) whose amplitude was 85% (95.5% in Hay et al., 2011) of the threshold current for inducing a somatic spike. In the middle of this stimulus, a 5-ms square pulse current was injected at the soma with an amplitude 40% (37.9% in Hay et al., 2011) of the threshold amplitude for eliciting a spike. Preceding or following the onset of this short pulse, an EPSP-like current was injected at the apical dendrite at a distance of 600 or 850 μm (700 μm in Hay et al., 2011) from the soma. The rise and decay time constants of this EPSP-like current were 0.5 and 5 ms, respectively, and its maximal amplitude was 50 or 15% (44.0% in Hay et al., 2011) of that needed for eliciting a somatic action potential. In the down state, the long square pulse current was absent, and to compensate for this, the somatic stimulus was suprathreshold with an amplitude 135% of the threshold amplitude for eliciting a spike (the exact value not given in Hay et al., 2011, but in Mäki-Marttunen et al., 2016 similar responses were obtained with an amplitude 136.4% of the threshold amplitude). The relative values were used in order to transfer the protocol from one model (the Hay model) to the others (the Almog model and the altered Hay and Almog models): the absolute values of the amplitudes are given in Table S4.

Figure 3 illustrates the temporal sensitivities of the L5PC models to combined somatic and apical stimuli, both when the neuron was in an up or down state (up state modeled simplistically, see above). These temporal sensitivity windows were different for different neuron models. In both up- and down-state protocols, coincident stimuli tended to produce a larger dendritic voltage response than non-coincident stimuli (Figures 3A,B). Consistent with the fact that the bursting models had smaller Ca²⁺ channel conductivities (which caused smaller SK currents) and (in Hay model) larger Na⁺ channel conductivities at the apical dendrite, the bursting models (Figures 3D–E, Hay-A₁, Almog) had a wider sensitivity window, i.e. they generated large responses even when the inputs were temporally more apart, than the non-bursting and mildly bursting models (Figures 3D–E, Hay, Almog-A₄). In Figures 3F,G, we quantified the sensitivity of the L5PC to coincidence of apical and perisomatic inputs by integrals of the sensitivity windows of Figures 3D–E for one variant (SCN1A variant for Hay models, CACNA1C variant for Almog models), and the corresponding data for all variants of Figure 2 are shown in Figure S5. Figure 3 shows that the variants affected both the sensitivity of the coincidence detection as well as the amplitudes of the dendritic response to the combined stimuli. As expected, the variant effects were typically more constrained in the neuron models that expressed little or no bursting behavior (Hay model, top panel of Figure 3F; Almog-A₅ model, bottom panel of Figure 3G) than in the other neuron models, where the temporal sensitivity windows were wider. The altered Hay models generally produced a burst of action potentials given the somatic stimulation alone (see Figure 1), which back-propagated to the dendrites and induced a Ca²⁺ spike. This is reflected as a lack of temporal window for the Hay-A₁ model in the down-state protocol (bottom panel of Figure 3D, dashed) and thus modest effects of the variants (bottom panel of Figure 3F, dim): a proper temporal window was also absent for both up- and down-state protocols in Hay-A₂ to Hay-A₆ models (data not shown). Extensive simulations of all variants of Table A2 (Supplementary Material) showed that the effects of the variants on the integrals of the temporal windows of the up states (Figures 3F,G, Figure S5) were (moderately) negatively correlated with the effects on the steady-state firing frequency of Figure S3 (correlation coefficient was between –0.3 and –0.7 in the unaltered and altered Hay and Almog models when considering all variants, and between –0.62 and –0.91 when considering only variants of Ca²⁺ channel and Ca²⁺ transporter-encoding genes; Table 3A). Similar results were obtained for the down states, where the correlations were weaker in Hay models but stronger in Almog models compared to the data of Table 3A (not shown).

Furthermore, we also studied the second proposed up- vs. down-state mechanism in particular by considering a scenario where the L5PC received noisy synaptic inputs, and the distinction between the up and down states was determined by the rate of inputs to the synapses. The modification of the excitatory synaptic drive is a more likely explanation for the difference in excitability between the up and down states than that of the inhibitory drive (Wilson, 2008). Therefore, we modeled the down state as a state where the frequency of inputs to the excitatory synapses was decreased by 30% (f = 0.504 Hz) from the average estimated value in vivo (see Hay and Segev, 2015), while in the up state the frequency of inputs to the excitatory synapses was increased by 10% (f = 0.792 Hz). These frequencies of synaptic activation led to firing rates 0.06 and 4.9 Hz in the down and up states, respectively, when no other stimuli were applied. The up-state firing frequency is in line with the experimentally observed values, as a range from 1.1 to 8.2 Hz with a mean of 3.5 Hz was observed in Neske and Connors (2016). For the down state, no systematic estimates have been reported, but the firing frequency has been observed to be extremely low compared to the up-state firing frequency, yet non-zero (Sanchez-Vives and McCormick, 2000; Lőrincz et al., 2015).

Figure 4 confirms the results of Figure 3 for control neuron and an SCN1A variant using an in vivo-like background synaptic firing and the description of up and down states as states with increased or decreased (respectively) rates of random synaptic inputs. The properties of the somatic current and the apical EPSP-like current were the same as in the experiment of Figure 3, except for the amplitude of the apical EPSP-like current in the up state, which was set 25% of the threshold amplitude to constrain the response. These simulations were performed only for the Hay model, and to show consistent trends, 125 repetitions of each simulation with different random number seeds (affecting the timing of presynaptic spikes) were run. Unlike in the noiseless Hay model but consistent with the Almog-model predictions, these results from the noisy Hay model show a wider temporal window in the up-state protocol than in the down-state protocol (Figure 4D). In both noisy and noiseless cases, the Hay model shows higher amplitudes in down-state protocols where the presynaptic spiking was sparser (noisy, Figure 4D) or the proximal dendritic depolarization was missing (noiseless, Figures 3A,D) but where the somatic stimulus was greater than in the up-state protocol. By contrast, in the Almog model and all altered Almog models, the down-state amplitudes were always lower than the up-state amplitudes (Figures 3B,E), despite an even larger compensation of the somatic stimulus amplitude (in the Hay model, the amplitude is 2.8 times larger in the down-state than in the up-state protocol, and in the Almog models, 3.4 times larger). These observations reflect the differences in the Ca²⁺-activated K⁺ currents between the two models, as discussed in Mäki-Marttunen et al. (2017). Importantly, the effects of the SCN1A variant were qualitatively the same in the noisy case (Figure 4) as in the noiseless case (Figures 3A,D)

Taken together, these data suggest that SCZ-associated variants can cause robust alterations of the communication between perisomatic and apical dendritic regions during both up- and down states of the L5PCs. These alterations include both amplification of the dendritic response to coincident stimuli at apical dendrite and the perisomatic region as well as widening of the temporal window of the inter-stimulus interval in which a large dendritic response is induced.

The medium/slow afterhyperpolarization current in the L5PC, mostly mediated by the Ca²⁺-activated K⁺ currents in L5PCs (Poirazi et al., 2003; Stocker, 2004), contributes to the L5PC adaptation and is thus an important regulator of the neuron firing behavior within a time scale of tens to hundreds of milliseconds following the spike (or other event that leads to large cytosolic Ca²⁺ transients). This time scale matches exactly the time scale of the PPI of the startle, which is impaired in SCZ patients (Turetsky et al., 2009). L5PCs are not a fundamental part of the mammalian startle network, but the insights gained from the biophysically detailed L5PC models can be used to make hypotheses on the cellular mechanisms in neurons within the startle network, such as the giant pontine reticular formation neurons (cf. Mäki-Marttunen et al., 2018a). In a similar manner as in Mäki-Marttunen et al. (2016), we studied how the variants of voltage-gated ion-channels and Ca²⁺ transporters affect the L5PC adaptation. To do this, we distributed 3,000 excitatory conductance-based alpha synapses (τ = 5 ms, E_rev = 0 mV) across the apical dendrite of the model L5PC (Figures 5A,B). We used the bisection method to find the threshold conductance g_th for making the neuron fire in response to simultaneous activation of these synapses. We then activated the synapses at time t = 0 with a suprathreshold stimulus (amplitude 10% above g_th) and again at time t = t_ISI with an amplitude c × g_th. Depending on the interstimulus interval t_ISI, the threshold factor c for eliciting an additional spike was smaller or larger than 1, reflecting a prepulse facilitatory or inhibitory effect, respectively. Figures 5C,D shows the threshold factors c as a function of the inter-stimulus interval ISI for different scalings of a CACNA1C variant according to the Hay model, Almog model, and altered Almog models. In Figures 5E,F, the maximal threshold factors c (across ISIs larger than 40 ms) are quantified for the CACNA1C variant of Figure 2, and Figure S6 shows the corresponding data for the rest of the variants of Figure 2.

Comparison of the variant effects on steady-state firing (Figure 2 and Figure S2) and prepulse-mediated adaptation (Figure 5 and Figure S6) indicates that the variants that increased the steady state firing decreased the adaptation. We confirmed this trend by simulating the effects of all variants of Table A2 (Supplementary Material) on these phenotypes. We found that the effects of the variants on maxima of the prepulse-inhibition curves of Figure 5 were (mostly) negatively correlated with the effects on the steady-state firing frequency of Figure S3 (correlation coefficient was –0.72 in the Hay model, 0.11 in Hay-A₁ model, and between –0.63 and –0.91 in the unaltered and altered Almog models; Table 3B).

As reported in Mäki-Marttunen et al. (2016), the Hay-model neurons express a relatively strong adaptation current, and the variants have a modest but observable effect on the amplitude of the prepulse-mediated adaptation measured as the firing threshold of a stimulus following a suprathreshold stimuls (Figure 5E and Figure S6A). This prepulse-mediated adaptation is mostly contributed by the SK current (cf. Mäki-Marttunen et al., 2018a), which is revealed by the comparison of the amplitude and kinetics of different potassium currents in the soma and the stimulated dendrite (Figures S7A,B): only SK current has the required strength and long decay time needed to cause such an inhibition. Here, we found qualitatively similar results for the Almog model neurons and altered Almog model neurons (Figure 5F and Figure S6B). In the Almog-model neuron, the dendritic SK currents are relatively large and long-lasting (Figure S7C), while the somatic SK currents are small (Figure S7D). In the altered Almog models, the magnitudes of the dendritic as well as somatic SK currents are increased (Figures S7E–J). This may explain why both the magnitude of the prepulse-mediated adaptation and the effects of the Ca²⁺ channel variants therein are larger in the altered Almog models than in the original Almog model (Figures 5E,F).

Taken together, our results suggest that the magnitude of an SK current-mediated adaptation by a prepulse stimulus can be crucially altered by SCZ-associated gene variants in L5PCs and that these predicted effects are qualitatively similar in many different L5PC models.

The analyses above shed light on how the neuron responses to stimuli of different types are altered by the SCZ-associated variants. Here, we aim at deeper characterization on how the SCZ-associated genes control the neural response to different combinations of spatially distributed stimuli. We divided the apical dendrite into six regions (0–200, 200–400, 400–600, 600–800, 800–1,000, and > 1,000 μm from the soma) and considered the basal dendrites as the seventh region. 1,000 excitatory conductance-based alpha synapses were uniformly distributed to each region, and as above, the threshold conductance for inducing a spike in the Hay model was determined (Table 2A). We then considered all possible binary combinations of input patterns, where each region is either non-stimulated (all synaptic conductances were set to zero) or simultaneously stimulated (all synaptic conductances were set to half of the threshold conductance). We determined the neuron output for each of these input combinations by considering seven different discrete measures: whether the neuron spiked (number of spikes induced; here 0–3) and whether the Ca²⁺ concentration in the middle of each of the six apical regions remained below 0.00011 mM (0), was between 0.00011 and 0.00012 mM (1), or above 0.00012 mM (2). These threshold values were chosen as approximate classifiers between low, medium, and high activation of intracellular Ca²⁺ signaling; depending on the outcome, the synapses in the region could be depressed (medium [Ca²⁺]) or potentiated (high [Ca²⁺]) as a consequence of synaptic inputs (Artola and Singer, 1993). Determining detailed, physiologically realistic values for these thresholds is out of the scope of this paper—instead, we here aim to show that the SCZ-associated variants may affect the neural coding both in terms of immediate outcome (number of spikes) and long-term, dendritic Ca²⁺-dependent effects (number of [Ca²⁺] thresholds crossed).

We first analyzed the input-output relationships in the control Hay-model neuron. The upper row of each panel in Figure 6A shows the discretized response of the control Hay-model neuron in terms of these seven discrete measures for all 128 combinations of stimuli. As expected, all output measures were positively (mostly weakly) correlated with the presence of the corresponding input (Table 2B). The number of output spikes was most correlated (correlation coefficient 0.43 in the control neuron) with the presence of basal stimuli, but it was moderately correlated (correlation coefficient 0.41 in the control neuron) also with the presence of the input stimuli in the apical tuft. The Ca²⁺ output at each apical dendritic region was typically most correlated with the presence of the stimulus at the same region, except for the most proximal apical dendritic region (0–200 μm), which was most correlated with the inputs in three most distal regions (correlation coefficient 0.36), and the fourth apical dendritic region (600–800 μm), which was most correlated with the presence of inputs in the neighboring third region (400–600 μm; correlation coefficient 0.53). The strong correlation between the [Ca²⁺] in the region 600–800 μm from the soma and the presence of inputs in the region 400–600 μm from the soma could be explained by the “hot zone” of Ca²⁺ channels which spans these two regions.

We next analyzed the effects of the SCZ-associated variants on these input-output relationships. The bottom row of each panel in Figure 6A shows the discretized response of the ϵ=12 CACNA1C variant of Figure 2, and the Hamming distances between the control and variant output patterns (i.e., the number of disagreeing data) are displayed at the left-hand side of each panel (average Hamming distance across all ϵ=12 variants shown in parentheses). The input-output correlations were qualitatively the same for the control neuron (Table 2B) and for the CACNA1C variant neuron (data not shown). Figure 6B shows the number of spikes (as in the upmost panel of Figure 6A) across all ϵ=12 variants. These data show that the variants had largest effect on integration of inputs in the basal dendrites and the second apical compartment (200–400 μm): Out of 32 input patterns where inputs were present in these regions, 11 caused substantial variation (i.e., more than five variants showed a differing output) in the output number of spikes. The corresponding numbers for other pairs of compartments were between 0 and 7. These numbers were, however, highly dependent on the differential scaling of the synaptic inputs at separate locations (see Table 2B): When we varied the scaling coefficient of the synaptic inputs from 0.1 to 1.0 (in Figures 6A,B, 0.5 was used), the relevance of different dendritic compartments varied as well (data not shown).

We next addressed the effects of the SCZ-associated variants on the output diversity of the Hay-model neuron. We defined the output diversity as the number of unique output patterns across the given input patterns. This is also a measure of the coding capacity of the neuron as it reflects the number of unique ways in which the neuron may respond to (and thus code) the inputs. As an example, although the control neuron in Figure 6A was given 128 different input patterns, it only produced 18 unique output patterns, as, e.g., the pattern where 2 spikes, high [Ca²⁺] in the 1st, 5th, and 6th apical compartments, a medium [Ca²⁺] in the 3rd compartment, and low [Ca²⁺] in the 2nd and 4th compartment were induced was obtained for 27 different input patterns. We calculated the output diversity for the variants of Figure 2 using absolute synaptic conductances ranging from 0.0025 to 0.5 nS—i.e., unlike in Figures 6A,B where synaptic conductance at each compartment was relative to the threshold conductance, here each compartment that received inputs had an equal synaptic conductance. This range of synaptic conductances spanned the discretized space of each output measure, meaning that for each output measure (number of spikes or [Ca²⁺] at a given compartment of the apical dendrite) there were combinations of inputs that produced all values of the measure for the control neuron. To smoothen the data, the output diversity was averaged over eleven repetitions, where a small depolarizing or hyperpolarizing somatic stimulus (whose amplitude ranged from –0.1 to 0.1 nA in steps of 0.02 nA) was applied.

Figure 6C shows the average output diversity as the function of the synaptic conductance. The curves show an expected shape, where the output diversity is close to zero for small input conductances (all combinations produce the same output, namely, no spikes and no significant elevations of [Ca²⁺]), increases for intermediate inputs, and decreases again for large input conductances (where many combinations produce a large number of spikes and large elevations of [Ca²⁺]). Nevertheless, the variants of Figure 2 alter the shape of this curve. The CACNA1C and ATP2A2 variants (ϵ=12 scaling) that increased the firing rate in Figure 2 also increased the output diversity for the large synaptic conductances in Figure 6C, and likewise, the SCN1A variant that decreased the firing in Figure 2 also decreased the output diversity in Figure 6C. However, the ATP2B2 variant increased the firing rate (Figure 2) but had little effect on output diversity (Figure 6C). We confirmed the relationship between the integral of the f-I curve and the average output diversity by calculating the correlation between these measures across all variants (Table 3C). We considered the output diversity below (synaptic conductances 2.5 to 32.5 pS) and beyond (synaptic conductances 35 to 50 pS) the approximate peak of the output diversity to separately take into account the neuron responsiveness to small and large-amplitude synaptic stimuli. The variant effects on f-I curve integral were strongly correlated with the effects on output diversity at large synaptic drive (correlation coefficient 0.73), but only mildly correlated with those on output diversity at small synaptic drive (correlation coefficient 0.31). Our earlier results showed that due to the strength of the SK currents, the variants of Ca²⁺-channel and transporter-encoding genes had qualitatively different effects than those of Na⁺ and HCN-channel-encoding genes (Mäki-Marttunen et al., 2017). We therefore also calculated the correlation coefficients separately for Ca²⁺-associated non-Ca²⁺-associated variants. These data showed that f-I curve integrals were moderately correlated with the output diversity at both small (correlation coefficient 0.64) and large (correlation coefficient 0.68) synaptic drives among the Na⁺, HCN, and K⁺-channel encoding variants. Taken together, these data suggest that variants that increase the excitability also increase the diversity of the neuron output: in other words, they make the neuron respond in a wider variety of ways (in terms of downstream spiking activity and modes of plasticity) to a given set of inputs.

We replicated the results of Table 2 and Figure 6 using the Almog model in Table S5 and Figure S8. The results are qualitatively the same, although the relative importance of different sections is different in the Almog model, and the neuron output was on average more sensitive to the effects of the model variants. In the same way as in the Hay model data, the intracellular Ca²⁺ transients were higher in the middle and distal parts of the apical dendrite than in the proximal apical dendrite, except for the most proximal part (0–200 μm), which had large Ca²⁺ transients in Hay model but extremely small in Almog model (Figure 6A and Figure S8A). In addition, both models predicted that the variations in the size of the Ca²⁺ transients caused by the SCZ-associated genes were on average larger in the middle of the apical dendrite than in the proximal and distal apical dendrites (both models predicted largest average d_H across the variants for the Ca²⁺ transients 400–600 μm from the soma, these values were 10.41 in the Hay model and 38.8 in the Almog model). In the Almog model, Ca²⁺ dynamics were described also in the basal dendrites, but the Ca²⁺ concentration never reached the smaller threshold in this region (Figure S8A). Unlike in the Hay model, the number of subsequent spikes (N_spikes>1) was highly sensitive to the effects of the variants in the Almog model (Figure S8B), causing a large average d_H (25.62 in the Almog model, 3.46 in the Hay model). Moreover, the correlation between f-I curve integral and predicted neuron output diversity across the variants was slightly larger for small input conductances (0.83) than for large input conductances (0.75) (Table 3C), which was opposite to Hay model predictions. Nevertheless, in a fashion similar to Hay model, the predicted output diversity had an inverted U-shaped curve with respect to the single-synapse conductance (Figure S8C). Except for the ATP2B2, HCN1, and KCNB1 variants, which had a minuscule effect on the output diversity in the Hay model (Figure 6C), the variants had qualitatively similar effects on the output diversity in the two models (Figure 6C and Figure S8C). These results confirm the possibility that the SCZ-associated variants may alter the neural coding capacity in L5PCs in terms of spiking and Ca²⁺-transient output.

Finally, we performed the simulations of the four considered phenotypes (steady-state firing, Figure 2; interaction of the perisomatic region and apical dendrite during up and down states, Figure 3; prepulse-mediated adaptation, Figure 5; and coding capacity, Figure 6) for combinations of variants to illustrate how small genetic effects could produce large effects when combined. From each gene, we chose the ϵ=12 variant that gave the maximal increase to the f-I curve integrals in the Hay model, except if all model variants of the underlying gene had negative effects (see Figure S3). Some of these model variants, namely those of CACNA1C, CACNB2, and CACNA1D, imposed cumulative effects on the same model parameters (see Table A2 in Supplementary Material), but most of them affected different types of transmembrane currents. Figures 7A–D shows the effects of this variant combination on the four phenotypes predicted with the Hay model, and Figures 7F–I shows the corresponding predictions from the Almog model. As expected, the effects of the variant combinations were larger than those of single model variants in Figures 2, 3, 5, 6. The coding by action potentials was also more dramatically affected by the combination of variants (Hamming distance between control and variant output patterns 16 (Hay model, Figure 7E) and 37 (Almog model, Figure 7J) than by the single variants (average Hamming distances 3.46 and 25.62; Figure 6A and Figure S8A). Importantly, the predictions of these phenotypes for the variant combinations followed the trends observed in Table 3: the combination of the variants, scaled as ϵ=12 (purple), increased the firing rate of the L5PC (especially for large inputs), suppressed the temporal window of sensitivity to paired apical and somatic stimuli, decreased prepulse-inhibiting ability of the L5PC, and increased its output diversity—da and vice versa for the ϵ=-12 scaling (green).

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