All models were run using the Topographica simulator 202122, which is freely available at https://github.com/ioam/topographica, with documentation at www.topographica.org. Topographica was developed for modeling the development of cortical maps, and has been typically used for computational modeling of aspects of low- and mid-level vision (e.g., orientation preference maps, orientation tuning, contrast sensitivity, aftereffects).

For this project, we used as our baseline the GCAL (gain control, adaptation, laterally connected) model 23, further described in 202122. GCAL incorporates features that are biologically realistic, including: 1) gain control at the lateral geniculate nucleus (LGN) level; 2) homeostatic adaptation of V1 responses based on a weighted sum of all inputs and limited by a logistic (sigmoid) nonlinearity; and 3) weights on excitatory and inhibitory lateral connections within V1, and on afferent connections to V1. These weights begin as isotropic (radially uniform), but subsequently modify in a self-organizing fashion upon repeated presentations of visual input and other forms of neural activity (e.g., lateral interactions) according to Hebbian (unsupervised activity-dependent) learning, with divisive normalization. Other characteristics of the GCAL model include: the assumption of single-compartment firing-rate neurons at the retinal ganglion cell, LGN, and V1 levels; hard-wired subcortical pathways to V1, including the main types of LGN neurons (e.g., On center-Off surround; Off center-On surround); roughly retinotopic projections from the retina to the LGN sheets to the V1 sheet; and separate parameters for excitatory and inhibitory activity (see Figure 1). GCAL has successfully modeled a wide range of phenomena expressed in a healthy V1, such as development of contrast-invariant orientation tuning and direction selectivity, development of ocular dominance, aftereffects, and surround suppression. Such modeling results show that these effects can be explained by a small number of canonical mechanisms 2124252627. These demonstrations also indicate that through visual experience, the statistical regularities of the environment are learned and encoded via the competitive processes inherent to Hebbian learning, and that these processes are constrained by gain control and homeostatic mechanisms to prevent runaway neural excitation and excessive plasticity in frequently activated circuits.

The basic GCAL model has four levels or sheets, each of which is implemented as a two-dimensional array of firing-rate neurons: a retina (24 x 24 density), LGN On and Off channel sheets (24 x 24 density), and a V1 sheet (48 x 48 density) (see Figure 1). Here, ‘density’ represents the number of simulated units (neurons) per unit area of the indicated sheet, corresponding to a square portion of the simulated visual field. The retina and LGN sheets thus have ¼ the number of units per visual area as does the V1 sheet. As illustrated in Figure 1, these four sheets are interconnected via sets of projections, or ‘connection fields,’ whose synaptic weights on the next sheet level are modifiable by Hebbian learning, and by changes in learning rates. Note that in GCAL models, the retina sheet is equivalent to the photoreceptor layer of the retina only. The model LGN sheet activity is an abstraction that represents all levels between the photoreceptors and the superficial layers of V1 that combine and transform the photoreceptor signals, including the retinal bipolar and ganglion cell layers and the LGN itself.

In all models presented here, activity was simulated in a series of time steps, with each corresponding to approximately 12.5 milliseconds in real time. At step 1 (t=0), the image reaches the retina, at step 2 the LGN On and Off sheets calculate their responses, at step 3 the output of the LGN On and Off sheets reaches the V1 sheet, and from steps 4 to 20 the activity within the V1 sheet propagates and settles through lateral inhibitory and excitatory connections. At the next step, the next stimulus is presented to the retina, and this repetition of sets of 20 steps continues for the number of iterations (i.e., stimulus presentation trials) chosen by the user.

For this project, we, for the most part, limited manipulations to those that provided the best fit to published data in our prior study 14, to determine if the same manipulations that led to CS and OR tuning changes could also account for increased neural variability. The specific manipulations we examined, and a rationale for their use, are described below.

Changes in Input Strength: Several lines of evidence indicate that afferent input to the LGN is reduced in schizophrenia. First, schizophrenia is associated with reduced retinal photoreceptor and bipolar cell signaling as indicated by smaller flash electroretinogram (ERG) waveform amplitudes 282930313233, and this change is more pronounced in patients at the outset of treatment for an acute psychotic episode than it is after several weeks of treatment. Preliminary data also indicate reduced activity in retinal ganglion cells 29, whose axons form the optic nerve and synapse directly onto the LGN. Second, reduced retinal signaling could be expected to reduce the strength of LGN output and therefore the strength of signals reaching V1, and indeed, reduced amplitudes of visual evoked potentials (VEPs) have been repeatedly observed in schizophrenia, with evidence that the problem is more pronounced in unmedicated patients (reviewed in 34,35,36). For this study, we examined V1 neural variability as a function of a 15% reduction in both retinal and LGN efferent signals. This manipulation led to compensatory hyperactivation in model V1 neurons in our prior modeling study 14, and provided a good fit to the increased CS seen in unmedicated first episode schizophrenia patients (reviewed in 3).

Reduced lateral inhibition within V1. The strength of inhibition within V1 was reduced based on findings of reduced GABA concentration in the visual cortex of people with schizophrenia 373839, as well as on findings of broadened OR tuning in cat after administration of the selective GABA_A antagonist gabazine 40. For the purposes of this study, a 10% reduction in V1 lateral inhibition strength was explored, as this change, in combination with an increase in the V1 afferent learning rate (see below, Changes in Plasticity) provided the best fit to published data on reduced CS and broadened OR tuning in chronic schizophrenia 14.

Two sources of evidence suggest elevated local excitatory activity within V1. One is the similarity between hyper-glutamatergic effects of ketamine administration in healthy volunteers and brain function in schizophrenia 414243. The second is evidence for elevated baseline gamma- and beta-band power and synchrony in people with schizophrenia 444546, suggesting abnormal network formation. For the purposes of this study, a 10% increase in V1 lateral excitation strength was explored, to parallel the extent of changes in inhibition.

Increased afferent learning rate at LGN-to-V1 connections. This manipulation to the Hebbian learning rate was based on findings of reduced reliability in neuronal co-activation patterns in schizophrenia 47, and tighter coupling between thalamic and cortical sensory processing regions in schizophrenia compared to healthy controls 48. This manipulation, in the form of a threefold increase in learning rate, was used only in combination with a 10% reduction in the V1 lateral inhibition rate, as this combination provided the best fit to published data on reduced CS and broadened OR tuning in chronic schizophrenia in our prior study 14.

We followed the model architecture and parameter settings from 23 except as specifically noted. Below, only the key algorithms related to activation, learning, and adaptation will be described, as these are most relevant to understanding the study’s schizophrenia-related parameter changes and findings relating to CS and OR tuning. These algorithms are the same as those described in 14.

In GCAL, the activation value (i.e., firing rate) in each unit in the retina sheet is the brightness of that pixel in the training image. Training stimuli (i.e., those that the models were exposed to during their 20,000 developmental iterations; see below Model Training and Development Strategy) were pairs of elongated Gaussians (see Figure 1). The center coordinates of each Gaussian were chosen from a random distribution, and the patterns covered an area that was calculated as 1/3 larger than the V1 sheet lengthwise and heightwise. This ensured that all regions of the V1 sheet were equally likely to be stimulated during training.

For each unit at position j in the LGN On and Off sheets (O), activation at time t +δt was calculated as:

Eq 1

Where γ_o = 14.0, set as a constant so that the overall strength of connections from the photoreceptor sheet to the LGN On and Off sheets will produce activations within the range of 0.0 to 1.0; γ_s is the strength of feedforward contrast gain control; Ψ_i is the activation of unit i in the two-dimensional array of neurons on the photoreceptor sheet from which LGN (On or Off) sheet unit j receives input (i.e., its afferent connection field F_j,P), and η_i,O(t) is the activation in other units in the LGN (On or Off) sheet at the previous time step (received over the suppressive connection field F_j,S); f is a half-wave rectifying function that ensures that activation values of LGN On and Off sheet units are always positive; the constant k = .11 ensures that output is always significant even for weak inputs; ω_ij is the fixed connection weight between retinal sheet unit i and LGN sheet j defined as a standard difference of Gaussian kernels (to ensure that an On sheet unit is activated only when there is both high center activity and low surround activity, and that an Off sheet unit is activated only when there is both low center activity and high surround activity) (23, see equation 3 for further details about derivation of these weight values); and ω_ij,srepresents the spatial profile of lateral inhibition received from other units in the same LGN sheet: these weights have a fixed, circular Gaussian profile (see 23, equation 4 for further details), and they have a divisive effect on LGN sheet activity that functions to implement contrast gain control.

For each unit (j) in the V1 sheet, activity is calculated as a function of the local connection field on and around j. The connection field (F_jp) consists of both incoming projections (p) from another (e.g., LGN) sheet (s_p), and activity in the inhibitory and excitatory lateral connections to that unit. The total contribution (C_jp) to j from projection p from each projection type is a dot product of the input and the weights in each connection field, or:

Eq 2

 Where η_i,p is the activation in unit i reflecting all of the neurons in V1 of a single projection type with which unit j is connected, and ω_ij,p is the connection weight from unit i to unit j in V1 for each projection type. The contributions from each projection type are combined to calculate the activation of each neuron j in V1 at time t as follows:

Eq 3

In the baseline GCAL model, the strength of each projection type was set as follows: γ_a(afferent) = 1.5, γ_E(excitatory) = 1.7, and γ_I(inhibitory) = -1.4. These values were shown to provide an appropriate balance between afferent and lateral influences, and between excitatory and inhibitory connections, and to lead to smooth and normally developing orientation and other topographic maps in the V1 sheet 23. The variable f is a half-wave rectifying function to ensure positive activation values, and its threshold value (θ) varies as a function of the average activity of the unit, as part of implementing homeostasis (see below, next two paragraphs).

From times t + 0.15 to t + 0.95 (i.e., steps 3 through 20 after each stimulus presentation), activity settles in V1, as noted above. At the end of this settling process, the level of activation in each V1 unit j is used to update the activation threshold point θ(see below) of each unit and to update the afferent and lateral inhibitory weights using the standard Hebbian learning algorithm (see below). V1 activity is then reset to 0.0 and a new stimulus is presented to the model.

The activation threshold θ determines how much a V1 unit will fire in response to an input. To determine the activation threshold for each V1 unit, a smoothed exponential average of its settled activity patterns is calculated:

Eq 4

where the parameter β, here set at 0.9999, determines the degree of smoothing in the calculation of the mean. The variable η̄_j is initially set to a target average for V1 unit activity, which is set to η̄_ja (0)= μ = 0.024. The threshold value θ is then updated as follows:

Eq 5

where λ is set at 0.01 and is the homeostatic learning rate, or the rate at which the average activity in each V1 unit is brought closer to the specified target value μ. If the mean activity level in a V1 unit begins to deviate from μ during model activity, the threshold θ is raised or lowered in order to bring the value back closer to μ.

As noted above, prior to any stimulus presentations to the model, connection weights from all projections to and within V1 (ω_ij,p) are roughly isotropic (radially uniform), with the constraint that weights for each type of projection are operative within different circular radii ((to enable, for example, a more narrow zone of lateral excitation surrounded by a larger zone where inhibitory effects operate, consistent with known cortical function49,50). Radii for GCAL, for afferent, lateral excitatory, and lateral inhibitory connections respectively were set to: r_A=0.27, r_E=0.1, and r_I=0.23.

As stimuli are presented to the model, afferent connection weights from the LGN On and Off sheets are updated each time the V1 settling process is completed (i.e., every 20 steps, or once per iteration). The updating is based on a Hebbian learning rule, leading to the development of connection strengths that represent correlations between the LGN On and Off unit outputs and the postsynaptic V1 response. At each iteration, the weight adjustment is therefore dependent on three factors: presynaptic (LGN) activity (η_i), postsynaptic (V1) activity (η_j), and the Hebbian learning rate (α), with specific influences as follows: 

Eq 6

where for each V1 unit j, α is the Hebbian learning rate for the afferent projection. Because unless Hebbian learning is constrained, it will lead to overly frequent changes in weight values and thus unstable network function 51, a constraint was added to impose divisive (postsynaptic) normalization on the V1 weight values. This well-understood mechanism 5253 was implemented via the denominator in the equation immediately above. For the baseline GCAL model, learning rates for afferent projections, lateral excitatory projections, and lateral inhibitory projections were set respectively to α_A=0.01, α_E=0.00, and α_I=0.03 as in23. As noted below, in several cases in this study modifications were made to these learning rates to explore the potential contributions of reduced stability within V1 on CS and OR tuning in schizophrenia.

For this study, we operationalized variability in three ways: 1) the spread of activation values (i.e., mean activation value of V1) across a large sequence of stimuli that varied only on a single dimension (i.e., variation due to stimulus change); 2) variability in activity across V1 neurons in response to a single stimulus (i.e., variation due to differential neuronal feature selectivity); and 3) fluctuation in the variability observed in case 2, across different levels of contrast. For practical reasons, we do not report on variability in response to a sequence of an unvarying stimulus because, due to characteristics of Topographica, even with biologically excessive modifications to models of V1 function, deviations from SDs of zero were uncommon and unreliable. On the other hand, reliable expansions or constrictions in firing rates were observed after selected manipulations when generating average responses across V1 to a sequence of 1000 stimuli that systematically increased in contrast (i.e., case 1 above), when examining responses to neurons within a model V1 in response to a single stimulus (i.e., case 2 above), or when examining changes in variability within V1 across the same basic stimulus at different levels of contrast (i.e., case 3 above). Because response variability is often observed in psychiatric populations, including schizophrenia 1054, across trials that differ in stimulus characteristics, this approach is relevant to much published data and to clinical observations.

We also recognize that in many ways, a high level of variability in firing rates across neurons within a single region of cortex is desirable in that it can reflect reduced redundancy. In the models we tested, however, we aimed to determine whether the specific manipulations we introduced would be associated with variability levels that were greater than normal, with normal being defined as the amount of variability previously demonstrated in published studies using the GCAL model in which multiple aspects of normal visual functioning were successfully simulated 14212324252627. We considered variability in excess of previously demonstrated normal levels to be pathological and to potentially be a cause of abnormal perception. This view is consistent with the ideas that precision in visual representation is the inverse of variability in sensory signals 55, and that excessive variability in brain activity in schizophrenia represents noise that is detrimental to efficient processing 56. For example, activation of orientation-selective cells whose orientation preference is far from that of a visible stimulus would lead to both increased variability within the distribution of V1 cells, and to broadened orientation tuning, which could then impact processes that depend on precise representation of orientation, such as contour integration, motion discrimination, and facial emotion recognition, all of which are impaired in schizophrenia 314. It has also been demonstrated that both decreased and increased dopamine D1 receptor stimulation causes less distinctiveness in representations of stimulus patterns, due to a higher than normal number of neurons being active in coding many patterns 7. Similarly, while some degree of neural noise can aid in detecting initially subthreshold signals by amplifying them so they can cross the detection threshold, an excess of noise can impair detection 68.

The basic model development strategy was the same as in 14. This involved: 1) developing a model of normal visual system function by training it with 10,000 pairs of stimuli (randomly-oriented Gaussians; see Figure 1); 2) making changes associated with a specific central nervous system (CNS) alteration (e.g., increased lateral excitation in V1; reduced retinal and LGN output, etc.) and running the model for another 10,000 iterations to allow for long-term adaptation to the new state; and 3) presenting new testing stimuli (sine wave gratings; see Figure 2) and measuring: a) average activation and variability across V1 neurons during a sequence of 1000 stimuli, b) variability within V1 neurons in response to a single stimulus, or c) the fluctuation in variability within V1 neuron firing rates in response to a single stimulus as contrast is increased (in separate models). In addition, to examine variability after shorter-term CNS changes, we also developed a normal visual system model to 10,000, implemented CNS changes, and then immediately (i.e., without 10,000 additional trials) presented new stimuli and measured activation and variability in the same three scenarios described above.

Post-development testing stimuli consisted of low or medium spatial frequency (SF) sine wave gratings (see Figure 2). For measurement of variability within a stimulus sequence, we used sequences of 1000 trials for each sine-wave stimulus, in which the contrast on each presentation increased by 0.001, starting at 0.000 and ending at 0.999, with contrast defined as the proportion of the possible input range 23. For each sequence of 1000 trials, data were collected on the mean and standard deviation (SD) of activity across all model V1 neurons. From these data, we calculated the relative standard deviation or coefficient of variation (CV) which is the standard deviation divided by the mean 5758. We used the CV as our primary measure of neuronal variability because while the SD typically increases as the mean increases [and so, increased variability can be secondary to increased activation in general [e.g., as in9,13the CV expresses variability independent of the mean, and so it is scale invariant. The CV can also be thought of as reflecting the volatility of observed values around an expected value. This metric is commonly used to characterize stock market investments, where it provides investors with an estimate of the risk associated with a given investment over short-term intervals. For measurement of variability between different feature selective neurons in response to a single stimulus, the mean, SD, CV, and excess kurtosis values of activation levels across V1 sheet neurons in response to a single sine wave stimulus were generated. As with calculation of variability across stimulus transformations, CV was used as our primary measure of variability, with excess kurtosis used to represent excessive activation in orientation selective cells whose preference is far from that of the presented stimulus.

Because increased excitatory activity can increase variability in neural network firing 5659, we expected that the greatest deviations from normal would arise from the manipulation that increased V1 activity to the greatest degree relative to unmodified models of V1 function. In our prior study this was a 15% reduction in retinal and LGN efferent activity, which led to a compensatory hyper-activation in V1 neurons. Because this effect was found mainly with low, but not medium SF stimuli, we expected to observe the same pattern here in terms of CV values. On the other hand, since impaired tuning is thought to be associated with increased variability as well, we expected that the model that best fit chronic schizophrenia data on reduced CS and broadened OR tuning in our past study (i.e., a 10% reduction in V1 lateral inhibition along with a threefold increase in the afferent learning rate at LGN-V1 connections 14) would also be associated with an increased CV, in addition to an increase in excess kurtosis values, as previously observed.

Table 1, Row 1 demonstrates that when a normal V1 model is developed over 10,000 iterations, and then a CNS change is made, followed immediately by presentation of a sequence of 1000 stimuli with incrementally increasing contrast: 1) the highest level of activation (over double that observed in the unmodified model shown in Row 2, Column 1), and the largest SD are observed in the model with increased V1 lateral excitation; while 2) the weakest level of activation but the highest CV values are observed in the model with decreased retinal and LGN output. In contrast, Row 2 demonstrates that when a normal model is developed over 10,000 iterations, and then a CNS change is made followed by 10,000 additional development trials in the new CNS state, followed by presentation of the same 1000 trial stimulus sequence: 1) the models with either increased V1 lateral excitation or decreased V1 lateral inhibition generated the lowest levels of activation (owing to the development of compensatory inhibition) but the largest CV value; and 2) the model with decreased retinal and LGN output generated the highest level of activation (owing to compensatory activity in V1, replicating Silverstein et al. (2017)), but the lowest CV value. Findings from these two models suggest that in the earliest phase of schizophrenia-related CNS dysregulation (which may begin in the prodromal phase), increased deviance from the normal range of excitation is likely to be related both to excess excitatory activity, as well as to variability that is independent of mean activity level (as reflected in CV values) that results from weakened input to V1. In contrast, with adaptation to these changes, a reversal occurs in which increased excitation is compensated for by increased lateral inhibition which lowers overall activation but leads to large increases in variability in V1 firing rates, while at the same time reduced input to V1 is compensated for by increased activation in V1 neurons but with reduced variability.

Single trial data are presented in Table 2. Row 1 demonstrates that when a normal V1 model is developed over 10,000 iterations, and then a CNS change is made, followed immediately by presentation of a single test stimulus: 1) the highest level of activation (over double that observed in the unmodified model shown in Row 2, Column 2), and the largest SD are observed in the model with increased V1 lateral excitation; while 2) the weakest level of activation but the highest CV and excess kurtosis values are observed in the model with decreased retinal and LGN output. In contrast, Row 3 demonstrates that when a normal model is developed over 10,000 iterations, and then a CNS change is made followed by 10,000 additional development trials in the new CNS state, followed by presentation of a single test stimulus: 1) the models with increased V1 lateral excitation or reduced V1 lateral inhibition generated the lowest levels of activation (owing to the development of inhibition) but the largest CV and excess kurtosis values; and 2) the model with decreased retinal and LGN output generated the highest level of activation (owing to compensatory activity in V1), but the lowest CV and excess kurtosis values (see Figure 3 and Figure 4). Taken together, and similar to the cross-stimuli data reported in the paragraph above, data from these two models suggest that in the earliest phase of acute schizophrenia-related CNS dysregulation, increased variability is likely to be related both to excess excitatory activity, as well as to variability that is independent of mean activity level (as reflected in CV values) that results from weakened input to V1. In contrast, with adaptation to these changes, a reversal occurs in which increased excitation is compensated for by increased lateral inhibition which lowers overall activation but leads to large increases in variability in V1 firing rates, while at the same time reduced input to V1 is compensated for by increased activation in V1 neurons but with reduced variability.

Data on CV values across a range of contrast levels (each multiplied by a factor of 2) from 20% to 80% contrast are shown in Figure 5, for the case where the model was run normally for 10,000 trials, followed by a model change and then immediately by presentation of a single LSF test stimulus at one of the contrast levels. It can be seen here that manipulations involving increased excitation or reduced inhibition generally follow the normal pattern of change in variability across contrast level. In fact, the overlap is so great among most conditions and the unmodified condition that most data series are occluded in the graph. However, the clear deviant condition is the one with reduced retinal and LGN signaling. No data point is plotted for 20% contrast in this condition because the mean was 0 and so the CV could not be defined. However, at 40% contrast (where there is excessive noise relative to the weak signal) the CV is excessively high, and it remains higher than in all other conditions even at 80% contrast, as noted above.

For the case where the model was run normally for 10,000 trials, followed by a model perturbation and then 10,000 additional trials in which the effects of the manipulation were increasingly manifested, see Figure 6. Here, where the system has adapted to reduced sensory input: 1) variability is lowest across all contrast levels in the condition of reduced retinal and LGN output, and variability is highest in conditions with either increased excitation or reduced inhibition.