Computational model links normalization to chemoarchitecture in the human visual system

A goal of cognitive neuroscience is to provide computational accounts of brain function. Canonical computations—mathematical operations used by the brain in many contexts—fulfill broad information–processing needs by varying their algorithmic parameters. A key question concerns the identification of biological substrates for these computations and their algorithms. Chemoarchitecture—the spatial distribution of neurotransmitter receptor densities—shapes brain function. Here, we propose that local variations in specific receptor densities implement algorithmic modulations of canonical computations. To test this hypothesis, we combine mathematical modeling of brain responses with chemoarchitecture data. We compare parameters of divisive normalization obtained from 7-tesla functional magnetic resonance imaging with receptor density maps obtained from positron emission tomography. We find evidence that serotonin and γ-aminobutyric acid receptor densities are the biological substrate for algorithmic modulations of divisive normalization in the human visual system. Our model links computational and biological levels of vision, explaining how canonical computations allow the brain to fulfill broad information–processing needs.


INTRODUCTION
Computational models aim to provide an explicit mathematical description of the information processing carried out by the nervous system.The visual system has long been used as a beachhead for these approaches (1).Computations first identified in the visual domain have later been observed in other sensory modalities or cognitive domains (2)(3)(4)(5)(6)(7)(8)(9).This has led to the proposal that some of these computations, such as receptive fields and divisive normalization (DN), may be "canonical, " i.e., that the same mathematical operations might be used by the brain in a variety of contexts (10,11).Depending on local parametrization of the algorithm (for example, the choice of weights in a receptive field), canonical computations can fulfill a variety of information-processing requirements and capture a wide range of seemingly disparate brain responses (11,12).A crucial question concerns the identification of biological mechanisms implementing and modulating the algorithms of these computations (13).
DN is considered a prime candidate for a canonical computation.Originally introduced to explain nonlinear properties of V1 neurons (2), it has since been applied in a variety of contexts (11).In vision, the population receptive field (pRF) is the region of visual space that elicits a response from a population of neurons (14,15).We have recently introduced a pRF model based on DN and shown that it unifies and outperforms existing pRF models by flexibly capturing a variety of brain responses (12).We have also shown that variations in specific model parameters explain different informationprocessing signatures observed throughout the visual hierarchy (12).
Chemoarchitecture-the spatial distribution of neurotransmitter receptor densities-is a major contributor to brain dynamics and functional properties in health and disease (16)(17)(18)(19).Multiple lines of evidence suggest the involvement of neurotransmitter systems in neural computations.Electrophysiological studies in animal models have shown that neurotransmitter receptors modulate neuronal responses in visual areas (20)(21)(22)(23)(24)(25).For instance, some receptors might influence stimulus sensitivity, by exerting differential modulations of response gain (e.g.5-HT1B and 5-HT2A); others might influence integration of visual stimuli over space, by modulating relevant neuronal properties, such as baseline activity and inhibition [e.g.5-HT1A and γaminobutyric acid (GABA)].Clinical studies in humans have also found altered neurotransmission and altered information processing in several brain disorders (15,(26)(27)(28)(29).However, explicit links between chemoarchitecture and brain computations have not yet been clearly identified.A key step in this sense is the inclusion of an algorithmic level of description; lacking this, any approach can identify, but does not necessarily explain, relationships between biological and functional levels (30).
Here, we address this issue by bringing together explicit mathematical modeling of brain responses (12,14) with chemoarchitecture data (31)(32)(33).Bridging the gap between biological and functional levels of description, we hypothesize that spatially varying receptor densities in the human cortex modulate the algorithms of canonical computations.This would explain how the brain can produce a wide range of responses and achieve a variety information-processing goals with a single computation.
On the basis of existing neurophysiological evidence (11,21,(23)(24)(25), we hypothesized the involvement of serotonin (5-HT1A, 5-HT1B, and 5-HT2A) and GABA receptors in DN.To test this hypothesis, we obtain algorithmic DN parameters from modeling of ultrahigh-field functional magnetic resonance imaging (fMRI) responses in a large cohort (n = 171) (12,14,34).We then compare DN computational model parameters with receptor density maps obtained from an independent positron emission tomography (PET) imaging experiment (31)(32)(33).This approach allows us to noninvasively investigate the algorithmic roles of receptor densities in cortical computations.We report a notable alignment between the fMRI-based computational model parameters and PET-based receptor density maps, consistent with our hypothesis.The pattern of association between the two independent datasets provides evidence in favor of the theoretical role of neurotransmitter systems as algorithmic modulators of canonical computations.In particular, we find strong support for the hypothesis that serotonin and GABA receptors modulate DN in the human brain.

DN pRF model links computational and biological levels of the human visual system
The properties of pRFs are robustly measurable with a variety of neuroimaging and recording methods (5,15,35,36).Distinct computational properties of pRFs, such as visual-field position, size, surround suppression, and compressive spatial summation (nonlinearity), have been observed, modeled, and shown to vary systematically throughout the human visual system (12,14,37,38).In this context, by surround suppression, we mean below-baseline deflections of time courses observed on the flanks of the central activation pRF, particularly prominent in early visual cortex (Figs.1E and 2A) (37); by compressive spatial summation, we mean sublinear scaling of responses to increasing input, particularly prominent in late visual cortex (Figs.1G and 2B) (38).We have previously shown that the DN model captures suppression and compression (as well as combinations of the two) via local variation in its algorithmic parameters, the activation, and normalization constants (Eq. 1 and Fig. 1B) (12).Here, following available evidence from neurophysiology (21,(23)(24)(25), we hypothesize that serotonin and GABA receptors might represent the biological substrate for algorithmic modulations of visual-spatial DN computations.In particular, we hypothesize that different densities of these receptors might underlie the observed variations in suppressive information-and compressive information-processing signatures, captured in the DN model by the activation and normalization constants (Fig. 1, D, E, and G) (12).Suppression, in the sense intended here, is likely to require inhibitory mechanisms, as well as an ongoing baseline activity that can be inhibited even in the absence of driving stimuli (39).Hence, mechanisms involved in the control of inhibition and baseline activity may be relevant biological correlates for suppression.Neurophysiological evidence suggests that GABA-mediated inhibition modulates responses to visual stimuli (24), while the 5-HT1A receptor has been recently shown to control baseline activity in rodent primary visual cortex (23).Hence, we hypothesize that different densities of GABA and 5-HT1A receptors might show opposite modulations of suppressive information-processing signatures (Fig. 1E).Following this hypothesis, we predict that the DN model activation constant should correlate positively with GABA receptor density and negatively with 5-HT1A receptor density in the visual system (Fig. 1F).
Compression implies relatively strong responses to weak input and relatively constant responses over a wide range of increasing input.Hence, mechanisms involved in the control of input and response gain may be relevant biological correlates for compression.Neurophysiological evidence suggests that 5-HT1B and 5-HT2A .the model describes responses to visual-spatial stimuli as the ratio of activation (input) and normalization (context).the modulation exerted by activation and normalization constants allows the dn model to capture a variety of information-processing signatures (12).At the implementation level, the dn model captures a combination of (C) circuitry, comprising feedforward, lateral, and feedback connectivity, and (D) receptors.(E) the activation constant of the dn model modulates the degree of suppression, shown in the difference between a pRF profile with strong suppression (dark blue line, high activation constant) and one with no suppression (green line, low activation constant).We hypothesize that inhibition and baseline activity may be relevant biological mechanisms for suppression.(F) Following this hypothesis, we predict that GABA and 5-ht1A receptor densities should have opposite correlations with estimates of the activation constant (positive and negative, respectively).(G) the normalization constant of the dn model (inversely) modulates the degree of compression, shown in the difference between an approximately linear response (black line, high normalization constant) and a strongly nonlinear (compressive) response (yellow line, low normalization constant).We hypothesize that modulations of input and response gain may be relevant biological mechanisms for compression.(H) Following this hypothesis, we predict that 5-ht1B and 5-ht2A receptor densities should have opposite correlations with estimates of the normalization constant (positive and negative, respectively).
receptors implement a bidirectional, complementary modulation of input and response gain (20,21).Activation of the 5-HT1B receptor has been found to boost neuronal responses to strong stimuli, while dampening responses to weak stimuli, whereas the converse pattern has been observed for the 5-HT2A receptor (20,21).Hence, we hypothesize that different densities of 5-HT1B and 5-HT2A receptors might show opposite modulations of compressive informationprocessing signatures (Fig. 1G).Following this hypothesis, we predict that the DN model normalization constant should correlate positively with 5-HT1B receptor density and negatively with 5-HT2A receptor density (Fig. 1H).

Population-level DN parameters are obtained from a large-scale fMRI dataset
To test our hypotheses, we obtained estimates of DN model parameters from the large-scale (n = 171) Human Connectome Project (HCP) retinotopy dataset (34,40).We then compare them with receptor density maps obtained by combining several independent PET imaging datasets (31)(32)(33).The DN model (12) recovers the well-known polar angle and eccentricity visual field maps throughout the human visual system (Fig. 3, A and C) (41).Fitting the DN model to each individual HCP participant data (Fig. 2) and averaging the resulting maps across participants, we also obtained groupaverage cortical maps of DN model activation and normalization constants (Fig. 3, B and D).We have previously shown that DN model constants systematically vary throughout the visual hierarchy, capturing the variation of suppression and compression across different visual areas (12).We replicate here our previous findings of systematic variation of DN model constants, in a larger cohort (fig.S1).For ease of visualization, we also show DN model parameter estimates in the visual system on flattened cortical surfaces (fig.S2).

Receptor densities vary systematically throughout the human visual system
Visual inspection of receptor density maps (Fig. 3, E to H) reveals that they covary with the functional visual hierarchy (Fig. 3, E to H).In particular, both 5-HT1B and GABA receptors show high density in V1, dropping starkly at the V1/V2 border; however, 5-HT1B and GABA receptors differ in their variation in later areas, as GABA density drops to its lowest values in dorsal areas, remaining relatively high in ventral visual areas; whereas the 5-HT1B receptor density drops to low values both in the dorsal and ventral directions.Conversely, the 5-HT1A receptor density is lowest in low-level visual areas and increases systematically in the posterior-anterior direction, reaching its highest values in temporal and ventral regions.Last, 5-HT2A receptor density displays a more complex pattern of variation in the visual system, with medium-high density in V1, decreasing in V2 and V3.5-HT2A receptor density then increases again, reaching highest density in temporal and lateral occipital regions, while instead decreasing to lowest density in late dorsal regions (Fig. 3, E to H).For ease of visualization, we also show receptor density maps in the visual system on flattened cortical surfaces (fig.S2).

Receptor densities correlate with DN model constants
Next, we computed correlations between the model constants and receptor densities (Fig. 4 and Materials and Methods) across predefined regions of interests (ROIs) of the HCP-multimodal parcellation (HCP-MMP) atlas (42), shown on flattened cortex in (fig.S2).
In sum, we find a notable alignment between the independently estimated receptor density maps (PET) and DN model constants (fMRI).In particular, we find direct evidence in favor of the hypothesis that 5-HT1A and GABA receptor densities represent biological correlates of the activation constant (suppression).We also find direct evidence in favor of the hypothesis that the 5-HT1B receptor  density might represent a biological correlate of the normalization constant (compression), while the role of the 5-HT2A receptor remains more ambiguous.

DN model constants exhibit multidimensional relationships with receptor densities
The one-dimensional correlations reported in Fig. 4 could be obscured by potential relationships between receptor distributions.To address this, we investigated whether the (linear) combination of receptor pairs that we hypothesized to be related to each DN model constant provides a better prediction compared to single receptors in the pair (Fig. 5).
We fit one-receptor and two-receptor general linear models (GLMs) on estimates of DN parameters obtained from one-half of the HCP participants and compare their predictions on parameter estimates obtained from the other, independent, half of participants, obtaining a cross-validated correlation coefficient wCC(cv).We find that the linear combination of GABA and 5-HT1A receptors provides a better cross-validated prediction for the activation constant [wCC(cv): 0.85] compared to either GABA alone [wCC(cv): 0.63] or 5-HT1A alone [wCC(cv): −0.67] (Fig. 5A).We also find that the linear combination of 5-HT1B and 5-HT2A receptors provides a better crossvalidated prediction for the normalization constant [wCC(cv): 0.58] compared to either 5-HT1B alone [wCC(cv): 0.44] or 5-HT2A alone [wCC(cv): −0.17] (Fig. 5B).Last, we compare the two-receptor models with surrogate models comprising, in turn, one receptor in the pair and randomized permutations of the other (Fisher's permutation test, 10 6 permutations) (Fig. 5, C and D).In all four cases, we find that the two-receptor models with both true receptor structures in the pair perform better than surrogate models with one true and one randomized receptor, in a statistically significant way (GABA, P < 10 −6 ; 5-HT1A, P < 10 −6 ; 5-HT1B, P = 2.70 × 10 −5 ; 5-HT2A, P = 4.74 × 10 −3 ).This result indicates that the improved performance of the two-receptor models is due to the combined spatial structure of the receptors in the pair and not simply obtainable by adding a regressor with an equivalent distribution of values.Notably, this result goes beyond the pairwise comparison results above by showing that the 5-HT2A receptor density also provides a plausible biological correlate for the normalization constant, consistent with our hypothesis, provided that the spatial distribution of the 5-HT1B receptor density is taken into account (Fig. 5D).
In sum, the multidimensional analysis consolidates the evidence for the hypothesized relationships between receptor densities and DN model constants.In particular, the combination of both receptors provides better cross-validated predictions for DN model constants, above and beyond individual receptors; statistical analysis confirms that the improvement is specific to the true spatial structure of receptor densities.Hence, these results provide additional evidence for the involvement of the hypothesized neurotransmitter systems in visual-spatial DN computations.Our findings indicate that GABA and 5-HT1A might represent independent and complementary correlates of the DN model activation constant and that 5-HT1B and 5-HT2A might represent independent and complementary biological correlates of the DN model normalization constants.

Principal components of receptor dataset correlate with DN model constants
Last, our hypothesis-driven approach may have overlooked other relationships between receptor density maps and DN model parameters.To examine this possibility in a data-driven way, we conducted a principal components analysis (PCA) of the receptor dataset and correlate these principal components with DN model parameters.For generality, we include all available receptors and all DN model parameters, including those for which we lack specific hypotheses.The full receptor dataset includes 5-HT1A, 5-HT1B, 5-HT2A, 5-HT4, 5-HTT, and GABA receptor density maps (31)(32)(33).We find that the receptor dataset is characterized by two major principal components, the key directions of variation in the receptor dataset.These two components together account for 81% of total variance and are the only components above equivariance threshold (Fig. 6A).The first component receives its largest contributions from GABA (negative), 5-HT1B (negative), and 5-HT1A (positive) receptor densities, with the 5-HT2A providing the smallest contribution.The second component receives same-sign contributions from all receptors, with the largest being the 5-HT2A receptor density and the smallest being the 5-HT1B receptor density (Fig. 6B).
Next, we computed correlations between the two major PCA components of the receptor dataset and each of the DN model parameters (Eq.1), i.e., activation and normalization pRF sizes (σ 1 and σ 2 ), eccentricity, polar angle, activation and normalization pRF amplitudes (a and c), and activation and normalization constants (b and d) (fig.S3).The first PCA component shows significant correlations with several DN model parameters; its strongest correlation is with the activation constant b (wCC: −0.83, P < 10 −6 ) (Fig. 6C).This result is notably consistent with our initial hypotheses: The first PCA component receives major contributions, opposite in sign, from GABA and 5-HT1A receptor densities and indeed correlates with the activation constant.The first component also receives a major negative contribution from the 5-HT1B receptor and indeed also significantly correlates negatively with the normalization constant (wCC: −0.62, P = 6 × 10 −6 ; fig.S3).The second PCA component is dominated by the 5-HT2A receptor, which is almost entirely absent from the main first component, suggesting a unique role of the 5-HT2A receptor in the visual system, orthogonal to the main component of variation.Notably, the normalization constant has the strongest correlation of all DN model parameters with this second PCA component (wCC: 0.34, P = 2.51 × 10 −2 ) (Fig. 6D).The direction of correlation is again consistent with the hypothesized relationship between 5-HT2A receptor and the normalization constant: The second PCA component receives a major (negative) contribution of 5-HT2A and indeed correlates with the normalization constant.Similarly to the relationship between the 5-HT2A receptor density and the normalization constant (Fig. 4D), this relationship appears to be bimodal; an association pattern is present in the combination of early, ventral, and temporal ROIs, while a cluster of dorsal ROIs follows a different pattern.
In sum, a PCA of receptor density maps provides additional evidence for the involvement of neurotransmitter systems in the suppressive and compressive modulation of DN computations throughout the visual system, in particular, for GABA and 5-HT1A as correlates of the activation constant and for 5-HT1B and 5-HT2A as correlates of the normalization constant.

DISCUSSION
By integrating explicit mathematical modeling of brain responses with chemoarchitecture data and drawing hypotheses from neurophysiological literature, here, we provide a new approach to bridging the gap between biological and computational levels of vision (30).We focus on DN, in the context of visual-spatial pRFs, to provide the link between these levels.In particular, following existing neurophysiological evidence, we hypothesized that the DN model's activation and normalization constants, which respectively capture suppressive information-and compressive information-processing signatures, might be explained by underlying variations in different receptor densities.To test these hypotheses, we compared estimates of DN computational model parameters obtained from a large-scale ) by principal components of the receptor dataset.We include in the PcA all receptors available in the dataset, even those for which we lack a specific hypothesis (5-ht1A, 5-ht1B, 5-ht2A, 5-ht4, 5-htt, and GABA).(B) composition of the first two PcA components of the receptor dataset, which together account for 81% of the variance.the first component receives its largest contributions from GABA, 5-ht1B, and 5-ht1A receptors, with the 5-ht2A providing the smallest contribution; conversely, the second component receives its largest contribution from the 5-ht2A receptor.(C) the first PcA component has its largest correlation with the activation constant, out of all dn model parameters.(D) the second PcA component has its largest correlation with the normalization constant, out of all dn model parameters.this data-driven analysis confirms the notion that neurotransmitter systems are involved in modulating visual-spatial computations and, in particular, that they might represent biological correlates of the dn model activation and normalization constants.
fMRI dataset, with receptor density maps obtained from an independent PET imaging experiment.We find evidence in favor of the hypotheses that 5-HT1A and GABA receptor densities are associated with the DN model's activation constant (suppression) and that 5-HT1B and 5-HT2A receptor densities are associated with the DN model's normalization constant (compression).These findings are supported by multiple analyses, thereby providing evidence of how varying receptor densities can allow the brain to fulfill a variety of information-processing requirements with a single canonical computation.Primary visual cortex (V1) is a major contributor to visual processing.V1 is an area serving a unique role in vision, both biologically and computationally, as the first cortical retinotopic area.On the biological side, V1 is the visual area with the largest surface area and has the highest/lowest densities for several receptors.In particular, the 5-HT1B receptor density drops sharply at the V1/V2 border (Fig. 3G).On the computational side, V1 has the highest DN model cvR 2 , which implies that its role is amplified in weighted correlations.It also has the highest or lowest values for several DN model parameters.We contend that these estimates should not be considered statistical outliers, i.e., an extreme value due to noise fluctuations of a uniform or normal distribution; rather, the extrema of DN parameters and receptor densities observed in V1 reflect its functional and biological uniqueness.The fact that independent datasets of DN model parameters and receptor densities both show extreme values in V1 provides additional support to the correspondence between DN computations and the underlying biology.To address the role of V1 in the relationships we identified, we repeated the analyses of Figs. 4 and 5 but in the absence of V1.We find that the results are near-identical, showing that the identified correlations are robust throughout the visual system, even in the absence of V1 (figs.S4 and S5).
Here, we compute correlations across predefined ROIs of the HCP-MMP atlas (42).We do this because (a) the effective resolution differs between the two datasets, complicating comparisons at the vertex level; (b) different ROIs are expected to be functionally distinct, and, hence, variation across ROIs should reflect global computational diversity, whereas we assume that similar computations exist within each ROI; (c) ROI-wide estimates provide more robust measures of both receptor densities and DN model constants than vertex-wise, as the latter suffer more the confounds due to warping and resampling, intrinsic in neuroimaging data; (d) for retinotopically defined ROIs such as V1, V2, or V3, within-ROI variance is by definition dominated by eccentricity, and, hence, vertex-wise correlations may come to include eccentricity-driven effects, unrelated from the hypotheses that we are testing here; (e) different spatial autocorrelations in the two types of data (DN parameters versus receptor densities) will affect statistical outcomes in ways that are difficult to fully control; and (f) consistency with previous studies of neurotransmitter receptors that also use ROI-wide estimates.Nevertheless, we acknowledge that DN model constants and receptor density maps may also vary within ROIs.To account for this, we also compute correlations on the vertex-wise data, with a dedicated approach to control for spatial autocorrelations (note S1).We find the results are near-identical, showing that the correlations are robust regardless of the choice of computing ROI-wide or vertex-wise approach (fig.S6).
Our primary results here are obtained by means of between-ROI correlations.In this context, within-ROI correlations may also of potential interest.However, investigating within-ROI correlations presents several limitations and confounds.In chiefly retinotopic ROIs such as V1, V2, and V3, within-ROI variance is primarily driven by eccentricity.Hence, any identified within-ROI correlation between receptor densities and other model parameters is likely confounded by eccentricity-driven variation.A theoretically viable solution to this confound could be to parcellate a retinotopic ROI by coarse eccentricity (e.g., foveal, parafoveal, and peripheral) and investigate correlations within each parcel.This approach is challenged by the limited stimulus extent of the HCP fMRI experiment and the relatively low spatial resolution and smoothing of PET data.The stimulus subtends 8° of visual angle, and, hence, we cannot obtain reliable estimates of DN model parameters beyond this eccentricity.This eccentricity range corresponds to a few centimeters of V1 surface; any subparcellation applied in this range quickly runs into the spatial resolution limits of PET data.In addition, the variation of DN model parameters within few centimeters of the same ROI is relatively small, since most of the relevant functional variation is between ROIs.Nonetheless, we note that our PCA (Fig. 6 and fig.S3) and our vertex-wise correlation analysis (fig.S6) do take into account within-ROI heterogeneity.The agreement of these analyses with our main results shows that the between-ROI analysis is robust with respect to within-ROI heterogeneity.Last, we note that in this work, we have taken a hypothesis-driven approach, informed by neurophysiological findings, to disentangle the wealth of potential associations between the two datasets.This approach does not allow us to exclude additional relationships between other pRF parameters and receptor densities.Eccentricity is certainly an interesting dimension in this sense.One noteworthy observation that we have not included in our main findings, since it was not part of our hypotheses, is a robust and statistically significant positive correlation between the 5-HTT (serotonin transporter) receptor density and visual eccentricity (fig.S8).This finding may be of potential functional relevance and worth of future investigation.
The association between the 5-HT2A receptor and the normalization constant is supported by the multidimensional and datadriven approaches.In particular, combining 5-HT1B and 5-HT2A receptor densities provides a better (cross-validated) prediction for the normalization constant compared to the 5-HT1B receptor alone, and the contribution of the 5-HT2A receptor is statistically significant compared to surrogate models comprising 5-HT1B and randomized permutations of the 5-HT2A receptor density (Fig. 5).In addition, the second PCA component of the receptor dataset receives the highest contribution from the 5-HT2A receptor and is found to have the highest correlation with the normalization constant out of all DN model parameters (Fig. 6).We note that the relationship between the 5-HT2A receptor density and the normalization constant appears to follow a bimodal distribution.The relationship appears to follow our hypothesis when considering early, temporal, and ventral areas of the visual system; a cluster of dorsal visual areas exhibits a quite different pattern of results.We speculate that the role of the 5-HT2A receptor in visual computations might differ across different visual streams.Specifically, the 5-HT2A receptor might play the role of the hypothesized visual-spatial normalization constant in lower, temporal, and ventral areas, whereas it might serve a different computational role in dorsal areas.We speculate that datasets involving richer stimuli, for example, probing temporal and multisensory aspects (43,44) of DN, may be needed to disentangle the computational role of the 5-HT2A receptor throughout the human visual hierarchy.
There are many mathematical formulations of DN that have been used in the literature (11,45).For example, some authors included additional parameters for exponentiation of activation and normalization components, while others do not.We argued that our choice of formulation strikes a fruitful balance between interpretability and parsimony (12).Nonetheless, we acknowledge that an equivalent model formulation can be obtained, without loss of generality, by fixing one of the four parameters a, b, c, and d to an arbitrary constant.For example, by setting d = 1 in the model Eq. 1.This formulation is advantageous in terms of parsimony, since it has one fewer parameter, but misses the straightforward interpretability of the normalization constant as modulator of compression.Estimates for the a, b, and c parameters in the d = 1 formulation are identical to those in our formulation, up to a divisive factor, i.e., the value of d at that cortical location.We find that the estimate of the activation constant b in this reduced d = 1 formulation (i.e., b/d in our formulation) still presents a positive correlation with GABA receptor density (wCC: 0.71, P < 10 −6 ) and a negative correlation with 5-HT1A receptor density (wCC: −0.59, P = 9.7 × 10 −5 ) (fig.S7).Correlations for the normalization constant cannot be computed, since the parameter is kept fixed to a constant value d = 1.This result proves that the proposed relationships between chemoarchitecture and DN model parameters are not exclusive to a single formulation and, for the remaining parameter (activation constant), can also be obtained in a more parsimonious alternative formulation.
The stimulus used in the HCP dataset does not vary in contrast, as is usually done to estimate DN models' neurophysiological contrast-response curves (11).Unlike these models, the DN pRF model aims to capture nonlinearities in visual space, e.g., compressive spatial summation (38), rather than in contrast.The standard moving-bar stimulus is shown at a variety of visual-field positions, orientations, and directions of motion (14).Hence, each cortical location is sampled with a variety of stimulus intensities due to its pRF position and size preferences.This process is akin, in space, to sampling at different contrast intensities for a contrast-response DN model.In a previous study (12), we varied bar sizes and speeds and found the same pattern of results in the DN parameters as with a single moving-bar stimulus.Whether the neural mechanism of the contrast-response nonlinearity is the same as spatial nonlinearity is an interesting question for future research but beyond the scope of this paper.
Our approach presents a number of methodological and conceptual limitations that should be emphasized.We are subject to the limitations of blood oxygen level-dependent (BOLD) fMRI.The BOLD signal does not allow distinguishing the activity of underlying neuronal subpopulations but only provides a proxy for overall activation, or deactivation, of the local neuronal population (~10 6 neurons), relatively (in our case) a stimulus-free baseline.Hence, here, we cannot distinguish the contributions of individual neurons or different cell types.Nonetheless, we argue that the scale of neuronal populations is an appropriate scale for description of computational processes such as spatial vision, which rely on the combined activity of large numbers of neurons.This approach is complementary to neurophysiological approaches, which provide high specificity and sensitivity in probing single-neuron responses at the expense of spatial coverage and overall response.We are subject to the limitations of PET.PET imaging on its own does not provide information as to how receptor densities assayed with radioactive ligands map onto underlying neuronal and receptor function nor to the particular cortical layer or cell type where they would be found.A strength of our approach is indeed the possibility of connecting PET-measured receptor densities to computations.As our interest here lies with the computations performed at the level of large neuronal populations, as measured by fMRI, we believe PET is the measurement method that is most apt for comparison with fMRI results.Alternative methods to measure receptor densities exist, such as postmortem mRNA expression.The receptor density maps we use here have been compared to mRNA measures and found the be in agreement throughout the cortex (31,32).Approaches focusing on modeling single-neuron responses might instead benefit from more fine-grained measurements of receptor densities, such as mRNA expression.We are limited by existing datasets to analyzing group-level relationships between visual computations and chemoarchitecture.However, pRF parameters and receptor densities are known to vary across individuals (31,46).Future work using PET and fMRI in a single cohort of participants might be able to address whether the same associations can be observed at the individual level.However, this approach is limited by the health implications of PET neuroimaging.A complementary direction aimed at answering these questions could involve the administration of neuropharmacological agents known to be active at specific receptors and measurement of resulting changes in pRF parameters.We only investigate a computational DN model of spatial vision.Hence, our approach is limited to the brain areas where, given the signal-to-noise ratio of the available fMRI data, we are able to robustly measure and model responses to visual-spatial stimuli (here, approximately one-sixth of the entire cortical surface).The relationships that we identify here need not generalize to different models, different sensory modalities, or other brain areas.Nonetheless, DN is considered a prime canonical computation candidate and has been identified in a variety of contexts, and pRF modeling approaches have been used in different sensory and cognitive modalities.Hence, our results beg the question of whether the relationships between biological systems and computational processes presented here are limited to spatial vision or apply more widely throughout the brain.We only investigate receptors that sit at the intersection of those with openly available PET data and those for which hypotheses can be formulated on the basis of existing neurophysiological evidence.Our approach is not exhaustive; other receptors beyond those investigated in this study are known to be involved in visual computations.Most notably, acetylcholine and N-methyl-d-aspartate receptors (47)(48)(49)(50).We cannot exclude that the information-processing signatures that we investigate here might involve additional or alternative mechanisms.Furthermore, it is certainly likely that the receptors that we investigate here are also involved in other kinds of computations (for example, on the basis of current results, we have speculated that the 2A receptor might subserve a different functional role in dorsal visual areas, compared to the rest of the visual system) and nonvisual processes throughout the brain.Our results also square well with studies on the receptor-nonspecific effects of serotonin in macaque visual cortex [e.g., (22)].The administration of serotonin results in subtractive and divisive change to visual response profiles, without affecting tuning preferences.These observation are thought to reflect a net effect between facilitation and suppression, mediated by different receptor classes all activated by serotonin.Our approach is inherently correlational.Our hypotheses are formulated on the basis of the findings of neurophysiological experiments in animal models, which generally involve direct, causal manipulation of neurotransmitter receptors via optogenetic or pharmacological interventions.Our findings hint at strong relationships between the results of computational modeling and biological systems and would be cemented by investigations with a causal approach, for example, by applying a direct pharmacological intervention in humans that would directly engage some of the neurotransmitter receptors that we find here to be associated with computational model parameters of DN.
In sum, in this work, we propose that spatially varying receptor densities contribute to the range of information-processing signatures observed throughout the human visual hierarchy and captured by the computational parameters of DN, activation, and normalization constants.We formulate hypotheses and predictions based on neurophysiological evidence and test them by combining explicit mathematical modeling of brain responses and chemoarchitecture data.Our results support the hypothesized roles of neurotransmitter receptors as algorithmic modulators of canonical computations in the human brain.In particular, we find that 5-HT1A and GABA receptor densities show opposite correlations with the DN model activation constant, indicative of suppressive information-processing signatures, while 5-HT1B and 5-HT2A receptors show opposite correlations with the normalization constant, indicative of compressive information-processing signatures.By providing an explicit link (the DN pRF model and its algorithmic parameters) between biological (receptor densities) and computational (suppressive information-and compressive information-processing) levels of description, our results provide insights into brain function and open routes for computational neuropharmacology in the human brain.

MATERIALS AND METHODS fMRI data acquisition and preprocessing
We obtain fMRI data from the HCP dataset, after it has been preprocessed according to HCP procedures and projected to a standard space (34).The HCP dataset comprises 7-T fMRI data for 171 participants.fMRI scans were collected, while participant viewed visual-spatial stimuli and performed an attention task at fixation.The stimulus conditions consisted of (a) a bar moving through the screen in eight different direction, (b) a clockwise rotating wedge, (c) a counterclockwise rotating wedge, (d) an expanding ring, and (e) a contracting ring.For each participant, two identical runs were collected for the moving-bar conditions; for all other conditions, only one run was collected.For additional details, see (34).Here, we use group-average time courses from all conditions to replicate the HCP-style analysis and show that the DN model recovers wellknown eccentricity and polar angle maps; we use individual-participant data, from the moving-bar conditions only, to construct groupaverage maps of DN model parameters, since (a) the presence of two identical runs for the moving-bar condition enables the calculation of cross-validated model performance, which cannot be carried out for the other stimulus conditions; (b) when not carrying out crossvalidation, averaging of at least two time courses markedly increases signal-to-noise ratio; and (c) wedges and rings are highly stereotyped stimulation pattern that do not allow probing pRF properties with the level of specificity enabled by the moving-bar stimulus, which, as it approaches receptive fields from a variety of directions, provides a more diverse set of responses and more robust estimation of pRF properties.

Model fitting
The DN model is defined as (12) In the above equation, the spatial dependence of Gaussian pRFs G ≡ G(x, y) and the spatiotemporal dependence of stimuli S ≡ S(x, y, t) are omitted for brevity, and we denote Hence, (x 0 ,y 0 ) is the response central location in the visual field, σ 1 and σ 2 are the sizes of activation and normalization pRFs, a and c are their amplitudes, and b and d are the activation and normalization constants.The DN model is fit following the same procedure as (12).Briefly, the original Gaussian pRF model is fit, with a grid and iterative stage, to obtain initial estimates of pRF size and position.These estimates are used as starting parameters for the DN model fit, which also comprises a grid and iterative fit stages.First, a grid stage is performed vertex-wise over the additional DN model parameters, while keeping the pRF size and position fixed to the Gaussian estimates; next, an iterative fit over all parameters is performed, including pRF size and position.At the iterative fit stage of both Gaussian and DN models, we fit a one-parameter hemodynamic response function based on the scanning probe microscope basis set.This allows adjusting for the potentially different speed of the hemodynamic response across areas and participant.Fitting the hemodynamic response function speed also ensures that our results are not due to different hemodynamic response properties but reflect distinct information-processing signatures of different areas.

HCP analysis replication
Following the same steps of the HCP analysis (34), we concatenate fMRI time courses for all available stimulus conditions (wedges, rings, and bars) and average data across participants.We carry out pRF model fitting on this population-averaged time courses to obtain estimates of classic visual field maps of polar angle and eccentricity.

Population-level cortical maps of DN model parameters
We first obtain single-participant maps of DN model parameters by fitting the DN model to the mean of the two available bar-pass stimulus runs, for each participant, at each cortical location.Averaging the two runs ensures high signal-to-noise ratio, required for robust parameter estimates.However, the variance explained (R 2 ) obtained from this fit is not cross-validated and therefore subject to potential overfitting issues.Hence, it cannot be taken as an unbiased measure of model performance.To address this, we compute the crossvalidated variance explained (cvR 2 ), as follows.We fit the model on each data run separately, evaluate performance on the other, independent run, and average the resulting two values.This gives us single-participant maps of cvR 2 , which provide an unbiased measure , G 1,2 ⋅ S ≡ ∑ x,y of model performance for each participant, at each cortical location.Next, we compute population-level parameter maps as the average of single-participant maps, weighted by the respective cvR 2 .For each participant, only cortical locations with a cvR 2 > 0 are included because negative values cannot be used as weights and because model fits with a cvR 2 < 0 cannot be considered reliable as they are likely the results of overfitting or highly noisy data.To obtain a population-level map of model fit quality, we compute the mean cvR 2 across the participants with cvR 2 > 0 at each cortical location.As these maps provide us with a measure of confidence in the parameter estimates, we use it for weighting in all subsequent analyses.Note that since a potentially different number of participants may have a cvR 2 > 0 at each cortical location, population-level parameter maps may receive contributions from a different number of participants at each cortical location.To take this into account, we only include, in subsequent analyses, those cortical locations where the map resulted from averaging of least a minimum number of participants, i.e., cortical locations where the number of participants with cvR 2 > 0 was higher than a minimum threshold value.To obtain this threshold in a data-driven way, we use Lins cross-entropy, a method often used in image processing for image segmentation/binarization. For our data, Lins cross-entropy gives a threshold value of 27.3 participants.Hence, we include in population-level parameter maps and subsequent analyses only cortical locations where at least 28 participants had a cvR 2 > 0.

Evaluating one-dimensional correlations between receptors and DN model constants
One-dimensional correlations are quantified using the (weighted) correlation coefficient, where the weight corresponds to the average cvR 2 in that ROI.This ensures that the importance of data points in the correlation coefficient computation is proportional to the confidence in parameter estimates.We assess statistical significance of each one-dimensional correlation by building a specific null distribution and performing a two-sided Fisher's permutation test.To do this, we obtain a null distribution for the independent variable by computing 10 6 randomized permutations of its true observed values; for each permutation, we compute the weighted correlation with the true dependent variable.The (two-sided) P value is then computed as the fraction of permutations in the null distribution that have a correlation greater in magnitude than the weighted correlation between the dependent variable and the true independent variable.This procedure provides the probability that the correlation between the true independent variable (here, receptor density) and the true dependent variable (here, DN model constant) could also have been obtained with a "surrogate" independent variable, with the same distribution of values as the true one but randomized spatial structure.

Evaluating two-receptor models
To compare two-receptor and one-receptor models, we use a crossvalidated, weighted correlation coefficient.We fit one-receptor and two-receptor GLMs on maps of DN model constants obtained from one-half of all participants and evaluate the models on the other, independent half.Cross-validation ensures that two-receptor model do not outperform one-receptor models purely due to overfitting and allows a fair comparison between models with different number of parameters.We assess statistical significance of partial correlations of each receptor in the two-receptor models by building a specific null distribution and performing a two-sided Fisher's permutation test.To do this, we compute 10 6 randomized permutations of the true observed values of each receptor in the pair, while maintaining the true values of the other receptor in the pair; for each permutation, we fit a two-dimensional GLM to the true values of the DN model constant and compute the weighted correlation coefficient.The (two-sided) P value is then computed, for each receptor in the pair, as the fraction of permutations in the null distribution that have a correlation greater in magnitude than the correlation between the DN model constant and the two-dimensional GLM with both true receptor densities in the pair.This procedures provides the probability that the correlation coefficient brought about by the tworeceptor models could also have been obtained by combining one of the two receptors with a surrogate receptor, with the same distribution of values as the other receptor in the pair but randomized spatial structure.

ROI definition
We use ROIs defined by the HCP-MMP atlas (42).

Supplementary Materials
This PDF file includes: Figs.S1 to S8 note S1

Fig. 1 .
Fig. 1.A model-based approach to visual computations and relationships between normalization and receptor densities.(A)the brain combines direct input and context to fulfill a variety of information-processing goals.dn is considered a prime candidate for this computation.(B) the algorithm of our dn pRF model (eq.1). the model describes responses to visual-spatial stimuli as the ratio of activation (input) and normalization (context).the modulation exerted by activation and normalization constants allows the dn model to capture a variety of information-processing signatures(12).At the implementation level, the dn model captures a combination of (C) circuitry, comprising feedforward, lateral, and feedback connectivity, and (D) receptors.(E) the activation constant of the dn model modulates the degree of suppression, shown in the difference between a pRF profile with strong suppression (dark blue line, high activation constant) and one with no suppression (green line, low activation constant).We hypothesize that inhibition and baseline activity may be relevant biological mechanisms for suppression.(F) Following this hypothesis, we predict that GABA and 5-ht1A receptor densities should have opposite correlations with estimates of the activation constant (positive and negative, respectively).(G) the normalization constant of the dn model (inversely) modulates the degree of compression, shown in the difference between an approximately linear response (black line, high normalization constant) and a strongly nonlinear (compressive) response (yellow line, low normalization constant).We hypothesize that modulations of input and response gain may be relevant biological mechanisms for compression.(H) Following this hypothesis, we predict that 5-ht1B and 5-ht2A receptor densities should have opposite correlations with estimates of the normalization constant (positive and negative, respectively).

Fig. 2 .
Fig. 2. DN pRF model captures distinct information-processing signatures at different cortical locations.(A) example data and model fits from an early visual cortex location (v1).the time course is characterized by low compression, visible in the similarly constant slope of the response profile to the bar passes, and high suppression, visible in the negative deflections of the blood oxygen level-dependent (BOld) time course on the sides of all positive responses to bar passes.the former is captured in the dn model by a high value of the normalization constant, and the latter is captured in the dn model by a high value of the activation constant.(B) example data and model fits from a late visual cortex location [Middle temporal/temporal Occipital (Mt/tO)].the time course is characterized by high compression, visible in the nonlinear response profiles to the bar passes, and low suppression, visible in the lack of negative, below-baseline deflections in the time course.the former signature is captured in the dn model by a low value of the normalization constant, and the latter is captured in the dn model by a low value of the activation constant.insets: cortical surfaces with asterisks indicating the respective location of the time courses (left), schematic of the stimulus sequence, comprising 8-bar passes in cardinal and diagonal directions (34) (top), two-dimensional pRF profile in the visual field, highlighting positive (red) and negative (blue) locations, model fit quality (R 2 ), and estimates of dn model constants (right).

Fig. 3 .Fig. 4 .
Fig. 3. Estimates of DN model parameters and receptor density maps from independent fMRI and PET datasets.cortical maps of classic retinotopic visual field parameters, (A) polar angle and (C) eccentricity, obtained by replicating the hcP fitting procedure(34). the dn pRF model recovers the known retinotopic features of visual field maps.cortical maps of the dn model modulatory parameters, (B) activation constant and (D) normalization constant, obtained by fitting the model to the data of each hcP participant individually (n = 171).the parameters robustly and systematically vary throughout the visual system (see also fig.S1).(E) GABA, (F) 5-ht1A, (G) 5-ht1B, and (H) 5-ht2A receptor density maps in the visual system, obtained from Pet imaging data(31)(32)(33).A first qualitative analysis already indicates that variations in receptor densities reflect known functional properties of the visual system. in particular, a marked drop in both GABA and 5-ht1B density is evident at the v1 border, while the two receptor differ in their variation at later stages of the different visual streams.the 5-ht1A receptor shows a systematic pattern of increasing density in the posterior-anterior direction, from early to late visual areas; while the 5-ht2A receptor shows a more complex pattern, with relatively high density in v1, decreasing in v2/v3, highest density in lateral-temporal areas, and lowest density in dorsal areas.Outlines of the first visual regions of interests (ROis) displayed on all surfaces are based on the hcP-multimodal parcellation (hcP-MMP) atlas (see also fig.S2)(42).

Fig. 5 .
Fig. 5. Receptor pairs provide bidirectional, complementary modulation of DN model constants.(A) the linear combination of GABA and 5-ht1A receptors provides a better cross-validated prediction for the activation constant compared to either GABA alone or 5-ht1A alone.this suggests that the distributions of both receptors combine to provide a bidirectional modulation of the activation constant.(B) the linear combination of 5-ht1B and 5-ht2A receptors provides a better cross-validated prediction for the normalization constant compared to either 5-ht1B alone or 5-ht2A alone.this suggests that the distributions of both receptors combine to provide a bidirectional modulation of the normalization constant.(C) in-plane views of the two-receptor GlM shown in (A), showing the variation of the activation constant as function of GABA and 5-ht1A receptors, after the other receptor in the pair has been taken into account.(D) in-plane views of the two-receptor GlM shown in (B), showing the variation of the activation constant as function of 5-ht1B and 5-ht2A receptors, after the other receptor in the pair has been taken into account.Statistical significance was assessed with two-sided Fisher's permutation test (10 6 permutations).

Fig. 6 .
Fig.6.PCA provides additional evidence in favor of hypothesized relations between receptor density maps and DN model constants.(A) variance explained (R 2 ) by principal components of the receptor dataset.We include in the PcA all receptors available in the dataset, even those for which we lack a specific hypothesis (5-ht1A, 5-ht1B, 5-ht2A, 5-ht4, 5-htt, and GABA).(B) composition of the first two PcA components of the receptor dataset, which together account for 81% of the variance.the first component receives its largest contributions from GABA, 5-ht1B, and 5-ht1A receptors, with the 5-ht2A providing the smallest contribution; conversely, the second component receives its largest contribution from the 5-ht2A receptor.(C) the first PcA component has its largest correlation with the activation constant, out of all dn model parameters.(D) the second PcA component has its largest correlation with the normalization constant, out of all dn model parameters.this data-driven analysis confirms the notion that neurotransmitter systems are involved in modulating visual-spatial computations and, in particular, that they might represent biological correlates of the dn model activation and normalization constants.