PLoS Computational Biology

PLoS Comput Biol. 13(11): e1005856

A hierarchical stochastic model for bistable perception

http://orcid.org/0000-0003-1918-7357Stefan Albert, Formal analysis, Methodology, Software, Visualization, Writing – original draft¹, http://orcid.org/0000-0002-3773-0440Katharina Schmack, Data curation, Funding acquisition, Investigation, Writing – review & editing², Philipp Sterzer, Conceptualization, Funding acquisition, Writing – review & editing², http://orcid.org/0000-0001-5791-6405Gaby Schneider, Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – original draft¹*

Wolfgang Einhäuser, Editor

1. Institute of Mathematics, Goethe University, Frankfurt (Main), Germany

2. Department of Psychiatry and Psychotherapy, Charité Universitätsmedizin Berlin, Germany

Technische Universitat Chemnitz, GERMANY

The authors have declared that no competing interests exist.

* E-mail: schneider@math.uni-frankfurt.de

© 2017 Albert et al

DOI: 10.1371/journal.pcbi.1005856

Published online: 20 November 2017

Abstract

Viewing of ambiguous stimuli can lead to bistable perception alternating between the possible percepts. During continuous presentation of ambiguous stimuli, percept changes occur as single events, whereas during intermittent presentation of ambiguous stimuli, percept changes occur at more or less regular intervals either as single events or bursts. Response patterns can be highly variable and have been reported to show systematic differences between patients with schizophrenia and healthy controls. Existing models of bistable perception often use detailed assumptions and large parameter sets which make parameter estimation challenging. Here we propose a parsimonious stochastic model that provides a link between empirical data analysis of the observed response patterns and detailed models of underlying neuronal processes. Firstly, we use a Hidden Markov Model (HMM) for the times between percept changes, which assumes one single state in continuous presentation and a stable and an unstable state in intermittent presentation. The HMM captures the observed differences between patients with schizophrenia and healthy controls, but remains descriptive. Therefore, we secondly propose a hierarchical Brownian model (HBM), which produces similar response patterns but also provides a relation to potential underlying mechanisms. The main idea is that neuronal activity is described as an activity difference between two competing neuronal populations reflected in Brownian motions with drift. This differential activity generates switching between the two conflicting percepts and between stable and unstable states with similar mechanisms on different neuronal levels. With only a small number of parameters, the HBM can be fitted closely to a high variety of response patterns and captures group differences between healthy controls and patients with schizophrenia. At the same time, it provides a link to mechanistic models of bistable perception, linking the group differences to potential underlying mechanisms.

Author summary

Patients suffering from schizophrenia show specific cognitive deficits. One way to study these cognitive phenomena works with the presentation of ambiguous stimuli. During viewing of an ambiguous stimulus, perception alters spontaneously between different percepts. Percept changes occur as single events during continuous presentation, whereas during intermittent presentation, percept changes occur at regular intervals either as single events or bursts. Here we investigate perceptual responses to continuous and intermittent stimulation in healthy control subjects and patients with schizophrenia. Interestingly, the response patterns can be highly variable but show systematic group differences. We propose a model that connects these perceptual responses to underlying neuronal processes. The model mainly describes the activity difference between competing neuronal populations on different neuronal levels. In a hierarchical manner, the differential activity generates switching between the conflicting percepts as well as between states of higher and lower perceptual stability. By fitting the model directly to empirical responses, a high variety of patterns can be reproduced, and group differences between healthy controls and patients with schizophrenia can be captured. This helps to link the observed group differences to potential neuronal mechanisms, suggesting that patients with schizophrenia tend to spend more time in neuronal states of lower perceptual stability.

Introduction

The phenomenon of bistable perception has fascinated researchers for a long time [1, 2, 3]. Recently, the description of response patterns to bistable stimuli such as the Necker Cube, Rubin’s vase or rotating spheres with switching rotation direction gained increasing interest in computational neuroscience [4, 5, 6, 7, 8]. By modeling dynamic changes of perception during viewing of one and the same stimulus, one aims at providing potential explanations for neuronal mechanisms underlying perception and perceptual changes and to identify related brain areas as well as potential dysfunctions, e.g. in schizophrenia [9, 10].

Interestingly, the response patterns to continuously shown bistable stimuli often share common properties [7, 11]. Typically, the distribution of intervals of constant perception, termed dominance times, is unimodal and right-skewed, and extremely short dominance times, i.e., rapidly fluctuating precepts, are rare [12, 13, 14]. The dominance times under continuous stimulation are therefore often modeled as Gamma distributed [15, 16, 17, 18, 19, 20]. The mean of dominance times can be highly variable across subjects [14, 20], whereas the coefficient of variation (CV) is often comparable [21].

In comparison to a continuous presentation, intermittent presentation of a bistable stimulus, i.e., by repetitive interruption of stimulation for short time periods, has been observed to stabilize the percept if the interruption period is long enough, typically longer than 0.7 seconds [15, 18, 22, 23, 24, 25]. In this case, dominance times get longer and can also show a certain degree of periodicity [14]. In addition, such stable phases with long dominance times during intermittent presentation can also interchange with unstable phases of rapid percept changes. Fig 1 shows examples of response patterns to continuous and intermittent presentation of a bistable stimulus from the dataset reported in [10].

Fig 1View larger versionFig 1. Examples of response patterns to a bistable stimulus.Response patterns to continuous (green, A,B) and intermittent (blue, C-F) presentation from the dataset reported in [10]. While the distribution of dominance times tends to be unimodal in the continuous case, stable and unstable phases interchange in intermittent stimulation. In addition, response patterns can be highly variable across subjects.

Modeling studies with elaborated mathematical models have been proposed that can explain a number of properties of bistable perception like the distribution of dominance times under continuous stimulation [8, 17, 19, 20, 21, 26] or cyclic behavior and the impact of the duration of the stimulus presentation on the dominance times in intermittent stimulation [14, 18]. One key ingredient of these models of bistable perception is typically a competition between neuronal populations that correspond to the different percepts [14, 17, 18, 20, 27]. In order to account for stabilized perception in intermittent viewing, the use of multiple timescales for memory traces of past perception has been proposed by [14] and [18].

Many such models require a high number of parameters in order to describe the variety of response patterns. As a consequence, they can often hardly be fitted to experimental data, especially in the typical cases when only a few dozen dominance times are observed. In addition, the majority of models focus either on continuous or on intermittent viewing. Interesting models that are applicable to both cases have been proposed by [14, 17, 18].

The relevance of a joint description of continuous and intermittent viewing is illustrated here on a dataset including responses of patients with schizophrenia and of healthy controls to continuous and intermittent presentation of a rotating sphere with ambiguous rotation direction reported earlier in [9, 10]. In [10], an enhanced alternation rate for the group of patients with schizophrenia during intermittent presentation was reported. Interestingly, when we analyzed previously unpublished data recorded in the same participants during continuous presentation, the opposite could be observed [Fig 2; the data was collected during an initial training run for which the experimental procedures but not the results are described in 10]. Due to the differences in patterns and time scales between continuous and intermittent presentation, the potential neuronal mechanisms underlying the transitions between the different response properties remain unclear.

Fig 2View larger versionFig 2. Alternation rates in control subjects and subjects with schizophrenia.During continuous presentation healthy controls (C) showed higher alternation rates compared to patients with schizophrenia (S) (left), while the opposite could be observed for intermittent presentation (right). Each grey dot indicates the perceptual alternation rate from one individual participant, colored diamonds and lines indicate group medians with 25%/75% quantiles. Two-sided Wilcoxon tests yielded p < .1 for both continuous and intermittent stimulation.

Therefore, we propose here a new model for the description of response patterns to bistable perception that links the observed behavior in continuous and intermittent stimulation to potential underlying neuronal processes. First, the model should be able to describe the high variety of both, continuous and intermittent stimulation within one model framework. Second, we use a minimal number of parameters in order to allow parameter estimation and model fitting to the typically short experimental data. This can then allow the statistical investigation of differences between clinical groups.

Note that strictly speaking, the term ‘dominance time’ refers to slightly different objects in continuous and intermittent viewing. While during continuous presentation, switches occur from a dominant to a suppressed percept (percept-switch), dominance times during intermittent presentation consist of multiple continuous presentation periods, and switches typically occur because of different perceptual choices at the onset of the presentation (percept-choice) [28]. In the present model, the observed sequences of dominance times are treated as conceptually similar. This simplification allows for a parsimonious model description in both continuous and intermittent viewing but may not fully capture the relation between the perceptual processes in the two regimes.

The remainder of the article is organized as follows. First, we use a simple Hidden Markov Model (HMM) that describes the observed perceptual processes with a few parameters. For continuous presentation, one state produces independent and identically distributed dominance times with a two-parametric distribution. For intermittent presentation, switching between stable and unstable phases requires two hidden states with short and long dominance times, respectively. The HMM has the advantage that it allows straightforward model fitting and data description with a minimal number of parameters. However, it remains descriptive and lacks relations to potential underlying mechanisms. Therefore, we link the HMM to a hypothetical underlying stochastic model. This model is termed here Hierarchical Brownian Model (HBM) and intends to describe aggregated underlying neuronal activity, producing the observed behavioral responses.

The HBM is based on two main ideas: First, it assumes that switching between percepts results from two conflicting neuronal populations [cmp., e.g., 18]. In order to minimize the number of parameters, this process is reduced to a simple Brownian motion with drift that fluctuates between two thresholds, where the first passage times indicate state changes [similar to 21]. For continuous presentation, one therefore requires only two parameters, i.e., the drift of the Brownian motion and the threshold. The distribution of the resulting first passage times—i.e., dominance times—is then the same as in the HMM, with a simple relation between the two HBM and the two HMM parameters. Second, in order to describe intermittent presentation in the same model framework, we use a hierarchical model. The idea is to describe the switching between stable and unstable phases that is typical for intermittent presentation by using an analogous threshold crossing mechanism of conflicting neuronal populations. Specifically, we assume a second pair of neuronal populations whose corresponding Brownian motion modulates the drift and threshold of the first population pair and thus causes switching between stable and unstable phases. We give a set of model assumptions under which the HBM parameters are comparable to the HMM parameters, thus allowing both model fitting to experimental data sets and potential relation to underlying mechanisms. The parameter estimation is straightforward using maximum likelihood and the HBM can reproduce both, the unimodal distribution in the continuous presentation and the bimodal distribution of dominance times in the intermittent presentation, including also various different response patterns. Moreover, it allows the identification of specific differences between the clinical groups in [10] and relates these to the hypothesized underlying processes.

Results

A simple Hidden Markov Model for perceptual responses

As a first step to describe the processes observed in bistable perception, we reduce data analysis to the dominance times, i.e., the times between reported changes of the percept. As described above, the distribution of dominance times tends to be unimodal in the continuous case, while stable and unstable phases interchange in intermittent stimulation.

We denote the dominance times by di, i = 1, 2, …, n. In the unimodal continuous case, we assume the di to be the realizations of independent and identically (i.i.d.) distributed random variables Di, i = 1, 2, …, n (Fig 3A), where the Gamma or the Inverse Gaussian (IG) distribution are suitable two-parametric distributions [17, 18, 19, 21]. For comparability with the HBM, we focus on the Inverse Gaussian distribution here. For the intermittent case, we assume a HMM with a stable and an unstable state, which are hidden and produce long and short dominance times, respectively (Fig 3B). This requires two parameters for the switches between states, and two parameters for the distribution of dominance times in each state.

Fig 3View larger versionFig 3. A simple HMM for bistable perception.(A) One state describes a unimodal distribution of dominance times under continuous presentation. (B) Two states (stable, S, and unstable, U) produce long and short dominance times under intermittent presentation.

Formally, let Y ≔ (Yi)i = 1,…,n describe a Markov chain on {S, U}, where S and U denote the stable and unstable state, respectively. Let pSU = 1 − pSS and pUS = 1 − pUU denote the transition probabilities. The dominance times (di)i∈{1,2,…,n} are assumed to be Inverse Gaussian distributed and conditionally independent given Y, with mean and standard deviation given by (μS, σS) for Yi = S and (μU, σU) for Yi = U.

Note that independence of dominance times is assumed here in continuous presentation. This assumption enables straightforward parameter estimation and is in agreement with the observation that serial correlations of dominance times are typically not reported [e.g., 29, 30]. However, weak long-term dependencies of dominance times reported under continuous presentation [31] cannot be reproduced in the HMM. As such long-term dependence was not observed in the majority of cases in the present data set, also showing no group differences, we use here the simple assumption of independence. In addition, the two used HMM parameters are sufficient to capture the main group difference in the response properties reflected in the alternation rate.

Parameter estimation and precision of parameter estimates

Continuous presentation. For parameter estimation in the continuous case, we simply estimate the parameters of the Inverse Gaussian distribution from the dominance times d₁, …, dn, assuming these are realizations of i.i.d. random variables. Let

denote the sample mean of the observed dominance times. For the Inverse Gaussian (IG) distribution [32] with mean μ and standard deviation σ, the maximum likelihood (ML) estimators are given by

[Image: pcbi.1005856.e002.jpg]μ^=μ¯dandσ^=μ¯d31n∑i=1n(1di-1μ¯d).(1)

The precision of these ML estimators is investigated by parametric bootstrap for the parameter values estimated from the 61 response patterns to continuous presentation in the sample dataset of [10]. For each parameter combination (μi, σi)i = 1,…,61 we simulated 1000 response patterns with length T = 240 as in the original data. We then compared the estimators [Image: pcbi.1005856.e003.jpg](μ^i,σ^i)i=1,…,61 with the parameter values underlying the simulation using the relative error (RE), defined e.g. for μ as [Image: pcbi.1005856.e004.jpg]RE(μ)≔|μ^-μ|/μ. The median relative errors for the 61 parameter constellations are shown in Fig 4A. Out of these, 54 (89%) showed estimation errors with median REs less than 0.25 (across the two parameters μ and σ, black). The remaining simulations (gray) showed only few percept changes, n < 20, as well as large CVs (σ/μ, Fig 4B).

Fig 4View larger versionFig 4. Precision of parameter estimates in the one-state HMM.For each of 61 parameter constellations in continuous presentation, 1000 simulations were performed with sample sizes as in the original data. (A) Median of the relative error (RE) for each parameter. (B) The CV as a function of the number n of simulated dominance times. Black points indicate constellations with mean RE across the parameters smaller than 0.25.

Intermittent presentation. In order to estimate the parameters in the two-state HMM in intermittent presentation, we use the Baum-Welch-Alghorithm [BWA, 33, 34, 35]. See section ‘Intermittent presentation: Baum-Welch-Algorithm’ in the Materials and methods for details on the BWA, the choice of starting values and specifications in the sample data set. For details on HMMs see, e.g., [36, 37].

In order to investigate the estimation precision of this approach, we again apply parametric bootstrap to the 61 parameter combinations estimated from the response patterns to intermittent presentation in the sample data set. Fig 5 shows the median errors obtained in 1000 simulations for every parameter constellation for T₁ = 1200 s and T₂ = 3600 s. For pSS and pUU the absolute errors (i.e., [Image: pcbi.1005856.e005.jpg]AE(pSS)≔|p^SS-pSS|) are presented due to the small values of the two parameters. For the time horizon of the data, T₁ (panel A), 50 of the 61 parameter combinations yielded average errors (i.e. mean median errors across all parameters) smaller than 0.25 (black). The remaining cases (gray) showed a large μU, i.e., less distinguishable stable and unstable distribution, or small sample size n ≈ 10. For the larger time horizon T₂ (panel B), almost all parameter combinations showed errors smaller than 0.25.

Fig 5View larger versionFig 5. Precision of parameter estimates in the two-state HMM.For each of 61 parameter constellations in intermittent presentation, 1000 simulations were performed. Log(median) of RE (for μS, σS, μU, σU) or of AE (for pSS, pUU) for T₁ = 1200 (A) and T₂ = 3600 (B). Black lines indicate constellations with mean errors across the parameters <0.25.

Application of the HMM to the sample data set

Here we use the described methods in order to apply the HMM to the sample dataset presented in [10] consisting of responses to continuous and intermittent stimulation obtained from each of 29 patients with schizophrenia and 32 healthy controls (for details on experimental procedures see Materials and methods). While response patterns were highly variable across subjects, the CV of dominance times (mean 0.79, SEM 0.04) was comparable to other studies reported in [21]. Serial correlation of adjacent dominance times of the same percept was typically small (mean of Kendall’s rank correlation [Image: pcbi.1005856.e006.jpg]τ¯=0.02), and statistically significant on the 5% level in less than 7% of the cases, which is about chance level. Concerning long-term dependence [cmp. 31], deviations from the assumption of independent dominance times were not observed in 81% of the cases, and no differences were observed between the experimental groups (p > .1, Wilcoxon test, for details on the analysis see Materials and methods). A correlation between the alternation rates in continuous and intermittent stimulation across subjects was not observed in either group, comparable to the results of [14].

Model fit. By fitting the HMM to response patterns in continuous and intermittent presentation as described above, the typical properties of the observed response patterns can be reproduced in simulations, including unimodal distributions for continuous presentation and changes between stable and unstable stages in intermittent presentation and a high variety of response patterns (Fig 6). For example, subject C shows rather regular stable phases, separated by unstable phases, while subject D shows an irregular response pattern, subject E shows only stable phases, while subject F shows almost only unstable phases. The parameter estimates of these example subjects are given in Tables 1 and 2. Note also that the response patterns of seven of the 61 subjects were described better by only one (stable or unstable) distribution than by the two-state HMM.

Fig 6View larger versionFig 6. Comparison of empirical response patterns to patterns simulated with the HMM.Examples of response patterns to continuous (green, panels A and B) and intermittent (blue, C-F) stimulation repeated from Fig 1 and corresponding simulations within the HMM (orange). The parameter estimates are given in Tables 1 and 2, respectively.

See full table Table 1. Estimated HMM parameters of the response pattens to continuous presentation shown in Fig 1A and 1B.

See full table Table 2. Estimated HMM parameters of the response pattens to intermittent presentation shown in Fig 1C–1F.

In addition to the good correspondence in the response patterns, no strong deviations could be observed from the model assumption of Inverse Gaussian distributed dominance times (Fig 7).

Fig 7View larger versionFig 7. Comparison of distribution of dominance times with theoretical distribution.The theoretical IG distribution in the HMM fitted to the empirical distribution of dominance times for continuous (A) and intermittent (B) presentation shown in Fig 6A and 6C.

Comparison of repeated trials. The HMM approach also allows studying the reproducibility of response parameters of subjects across multiple sessions. To this end, we used an additional dataset from 105 healthy individuals that contained two separate sessions of continuous presentation [see 9, here we used the first two training runs from Behavioral Experiment 2 as described in this previous work]. These showed highly reproducible response patterns, i.e., a high correlation of parameter estimates of the IG distribution of the same individuals across different sessions (Fig 8).

Fig 8View larger versionFig 8. Reproducibility of response patterns.Parameter estimates of log(μ) (A) and log(σ) (B) of the IG distribution in two sessions with the same individuals, dataset reported in [9]. The logarithm was applied due to asymmetric distributions of the parameter estimates. Stars indicate highly significant (p < .0001) correlation of parameter estimates across different sessions.

In addition, we also used a likelihood ratio test [38] to investigate the null hypothesis of equality of the parameters of two Inverse Gaussian distributed samples with sample sizes n₁ and n₂, i.e. H₀: μ₁ = μ₂ and σ₁ = σ₂. The likelihood ratio derives as

with [Image: pcbi.1005856.e016.jpg]Si=∑j=1ni(dij-1-μ^i-1) for i = 1, 2, [Image: pcbi.1005856.e017.jpg]S3=n1/μ^1+n2/μ^2-n2(n1μ^1+n2μ^2)-1, n = n₁ + n₂ and S = S₁ + S₂ + S₃. Under H₀ the quantity [Image: pcbi.1005856.e018.jpg]Qn*≔-2(1-1/6[1/n1+1/n2]-1/[12n])logQn is approximately chi-square distributed with two degress of freedom. Thus, the test rejects the null hypothesis at level 5% if [Image: pcbi.1005856.e019.jpg]Qn* exceeds the 95%-th-quantile of the χ²(2)-distribution. In the sample data set, the likelihood ratio test did not reject the null hypothesis of equal parameter sets in 83 out of 105 subjects (about 79%). For a comparison, we performed 10000 permutations by randomly assigning a first trial of one subject to a second trial of another subject and performing the likelihood ratio tests on the permuted data sets. In the mean, the null hypothesis was not rejected in only about 36% of the randomly assigned pairs, with a maximum percentage across all permutations of 51%.

In summary, the response patterns of the same subject across multiple sessions showed a high degree of reproducibility, with a Pearson correlation coefficient of up to r = 0.76 between log(μ₁) and log(μ₂) (Fig 8). The similarity of response patterns for the same subject across multiple sessions was significantly higher than the similarity of response patterns between subjects.

Group differences. Finally, the HMM provides a relation between the underlying model parameters and the observed group differences reported in the introduction and in [10]. In the continuous case, the decreased alternation rate in the patients with schizophrenia is simply reflected in an increased mean dominance time [Image: pcbi.1005856.e020.jpg]μ^ (Fig 9A) in the one-state HMM. For intermittent presentation, we had observed an increased alternation rate in the patients with schizophrenia. The interpretation of this observation was not obvious due to the high variability of response patterns and particularly due to the fluctuation between stable and unstable state. The HMM provides a first explanation of this phenomenon by capturing important response properties in the parameter estimates, which showed the following group differences: In particular, the expected relative time spent in the stable state,

was higher in the control group (Fig 9B, detailed formula given in Eq (13) in the Materials and methods). As the main variable contributing to this difference, we observe that the probability [Image: pcbi.1005856.e022.jpg]p^SS to stay in the stable state was higher in healthy controls. In addition, the mean dominance time [Image: pcbi.1005856.e023.jpg]μ^U in the unstable state was slightly larger in the patients with schizophrenia. The degree of statistical significance was highly similar to the one reported in [10] (p < .1 for [Image: pcbi.1005856.e024.jpg]μ^, [Image: pcbi.1005856.e025.jpg]μ^U, [Image: pcbi.1005856.e026.jpg]p^SS and [Image: pcbi.1005856.e027.jpg]φ^S, two-sided Wilcoxon test). In the following section, we present a model that links these observed differences to potential underlying neuronal mechanisms.

Fig 9View larger versionFig 9. Differences in the HMM parameter estimates between subjects with schizophrenia and control subjects.(A) [Image: pcbi.1005856.e028]μ^ during continuous presentation and [Image: pcbi.1005856.e029]p^SS and [Image: pcbi.1005856.e030]φ^S for intermittent presentation (B). Each grey dot indicates one individual participant’s data, colored diamonds indicate medians, horizontal bars indicate 25%/75%-quantiles.

A hierarchical Brownian motion model

As described in the previous section, the HMM captures a high variety of response patterns both in continuous and intermittent viewing, including uni- and bimodal distributions of dominance times with alternations between stable and unstable states and a high variability across subjects. With its small number of parameters, the HMM can be fitted also to short data sections available empirically and therefore also capture differences between experimental groups.

However, the HMM description remains phenomenological and does not provide insight into potential neuronal processes. Also, it cannot provide explanations for potential effects that different lengths of blank displays could have on the response patterns, as discussed for example by [14, 22, 24, 25]. In addition, the HMM cannot represent the following interesting empirical observation: Before changing from stable to unstable state, the last dominance time tends to be shorter. Therefore, we introduce here a new model, called Hierarchical Brownian Model (HBM), which provides a potential link between the phenomonological description of the response and potential underlying neuronal processes. The HBM assumptions can also provide hypotheses on the effects of different lengths of blank displays and naturally yields shorter dominance times before a state change to the unstable state.

The HBM assumes two competing neuronal populations which indicate perception of right and left rotation, respectively. As has been proposed by various authors [14, 18], we implicitly assume mechanisms of self-excitation, cross-inhibition and adaptation across these neuronal populations, without explicitly modeling them in order to reduce the number of parameters and to allow for model fitting to short trials. In order to obtain a parsimonious model description, we again assume independence of dominance times by neglecting potential mechanisms of week long-term adaptation [31]. For possible model extensions compare section ‘Applicability and model extensions’ in the discussion. We use the simplified assumption that perception arises from the difference in the activity of the two populations, which is modeled here by a Brownian motion with drift [similar to 21] that fluctuates between two thresholds, where the first passage times indicate state changes. This results in two parameters for the case of continuous presentation that are directly linked to the two parametric distribution of dominance times in the HMM. Further, we describe switching between stable and unstable states in intermittent presentation by applying an analogous mechanism, which leads to a hierarchical model. We assume another hierarchical layer of neuronal populations and a corresponding Brownian motion which modulates the drift of the first population pair and thus causes switching between stable and unstable phases.

Continuous presentation

The HBM in continuous presentation (HBMc) simply assumes a Brownian motion with drift ±ν₀ between two borders, ±b, where the first hitting times of the borders indicate a percept change and lead to a sign change in the drift. As a potential neurophysiological interpretation, b could be considered the size of the activity difference between the L and R population required for a perception change and thereby be related to the respective population sizes. Roughly speaking, the speed of the drift ν₀ could be considered related to the inverse of the connection strengths within and across populations that engage in self excitation and cross inhibition.

Formally, let b > 0 be a fixed border, ν₀ > 0 be a drift and T > 0 a time horizon, and let (Wt)t∈[0,T] be a standard Brownian motion. The perception process P ≔ (Pt)t∈[0,T] is then defined by

[Image: pcbi.1005856.e031.jpg]dPt = Stν₀dt + dWt,  where P₀ =  - b, 

and the process St ≔ S(Pt, t) takes the value −1 if Pt last hit b and 1 if Pt last hit −b, with S₀ ≔ 1. For t ∈ (0, T] let

[Image: pcbi.1005856.e032.jpg]t*  ≔  t*(t)  ≔  sup{x|x < t, |Px| = b},  where t*(0) = 0

be the last time before t that Pt hit either b or −b. Then

[Image: pcbi.1005856.e033.jpg]St  ≔  S(Pt, t)  ≔  - sgn(Pt*).

The perception at time t ≥ 0 takes the value L (left) if St = 1 and R (right) if St = −1 and switches at the first-hitting times (Hi)i of the borders ±b defined by H₀ ≔ 0 and

[Image: pcbi.1005856.e034.jpg]Hi  ≔  inf{t ≤ T|t > Hi-1, Pt = Stb} i = 1, 2, …(2)

An example of such a process is shown in Fig 10. Panel A shows the process P, where the sign of the drift changes at each first hitting time of b or −b indicated by the process (Hi)i, which also marks switches in the percept (panel B).

Fig 10View larger versionFig 10. The HBMc.The first passage times (Hi) (green) of a Brownian motion Pt (black, panel A) with drift ±ν₀ at borders ±b indicate the times of the percept changes (orange, panel B). The Brownian motion is assumed to summarize the activity difference of two conflicting neuronal populations with only two parameters.

Parameter estimation in the HBMc makes use of the fact that the resulting dominance times

[Image: pcbi.1005856.e035.jpg]di  ≔  Hi - Hi-1,  i = 1, 2, …(3)

are independent and IG distributed, with a known relation between the HBMc parameters b and ν₀ and the parameters μ and σ of the IG distribution, given by μ = 2b/ν₀ and [Image: pcbi.1005856.e036.jpg]σ=2b/ν03 [32]. Note that the CV of the dominance time distribution is therefore given by 1/(2bν₀). Thus, an increase in the border b also increases the mean dominance time and decreases the CV. An increase in the drift ν₀ decreases the mean dominance time, while again decreasing the CV. The ML estimators [Image: pcbi.1005856.e037.jpg]b^ and [Image: pcbi.1005856.e038.jpg]ν^0 can be derived via

where [Image: pcbi.1005856.e040.jpg]μ^ and [Image: pcbi.1005856.e041.jpg]σ^ are derived from Eq (1). Explicit expressions are given in the Materials and methods.

The precision of parameter estimates is directly comparable to the HMM for continuous presentation.

Intermittent presentation

In the Hierarchical Brownian Motion model for intermittent presentation (HBMi), we require mechanisms for long dominance times in the stable state as well as for short dominance times in the unstable state. In order to describe the responses to intermittent and continuous presentation in one model framework, we assume the identical perceptual process as in the HBMc during phases of stimulus presentation. The periods of blank display represent the only difference to continuous presentation. In these periods, we assume additional neuronal mechanisms. In particular, we assume that the perceptual process then takes on one of two mean drifts, νS in the stable state and νU ≥ νS in the unstable state, with potentially opposite signs of ν₀ and νS for increased stability (Fig 11). Note that the drifts νS and νU are not necessarily constant across the whole period of blank display, but they denote the mean drift of the process, which is sufficient to describe the distribution of dominance times. Interestingly, additional assumptions on the temporal behavior of the drift terms could also allow describing the impact of the lengths of blank displays (cmp. Discussion). Further, in the unstable state the border bU at which perception and drift direction change is assumed smaller than the border bS during stable perception. Switches between the stable and unstable state will be caused by a similar mechanism in a so-called background process B described later in this section.

Fig 11View larger versionFig 11. The perception process P in the HBMi during intermittent presentation.During presentation, P has drift ν₀. During blank displays (yellow), P has drift νS in the stable phase (A), and drift νU in the unstable phase (B). Typically, we have νS ≤ ν₀ and νU ≥ ν₀. The borders are bS (light blue horizontal line) in the stable phase and bU (blue horizontal line) in the unstable phase.

Within a state (S or U), the fluctuation of the perception process between the borders is assumed analogous to the HBMc, except that the borders are dependent on the hidden state and that the drift is ν₀ during presentation and νS or νU during blank display. Formally, we denote by PR and BL the sets of all periods of stimulus presentation and blank display, respectively. Assuming that we start a trial with a presentation interval and then switch regularly between presentation intervals of length lp and blank display of length lb, PR and BL are given by

The perception process (Pt)t is then given by

where [Image: pcbi.1005856.e044.jpg]Y˜t∈{S,U} denotes the hidden state and (Wt)t denotes a standard Brownian motion. As a result, the mean drift per second is given by

[Image: pcbi.1005856.e045.jpg]νS*≔lb·νS+lp·ν0lb+lpandνU*≔lb·νU+lp·ν0lb+lp(4)

for states S and U, respectively. Because the periods lb and lp are typically short in relation to a dominance time, the behavior of P can be approximated by a Brownian motion with absolute drifts [Image: pcbi.1005856.e046.jpg]νS* and [Image: pcbi.1005856.e047.jpg]νU*, respectively. As in the HBMc, the sign of the drift St ≔ S(Pt, t) changes at every first hitting time of the respective border, i.e.,

We initialize [Image: pcbi.1005856.e049.jpg]P0=-bY˜0 for the initial state [Image: pcbi.1005856.e050.jpg]Y˜0 which is the stable state with probability [Image: pcbi.1005856.e051.jpg]πS≔P(Y˜0=S). The perception then takes the value L if St = 1 and R if St = −1 and switches at the first-hitting times (Hi)i of the borders ±bi comparable to Eq (2). Note that perception also changes during blank display. The dominance times are therefore again given by di ≔ Hi − Hi−1, i = 1, 2, ….

In order to describe the switching between the two states S and U, we use an analogous upper hierarchical level with another pair of conflicting neuronal populations. Their difference activity is described by a so-called background process B ≔ (Bt)t (Fig 12A, middle panel). B is also assumed to be a Brownian motion with drift. Its drift is assumed to vanish during presentation and to take the value νB > 0 during blank display, i.e.,

[Image: pcbi.1005856.e052.jpg]dBt={dW˜t,ift∈PRνbdt+dW˜t,ift∈BL,Y˜t=S-νbdt+dW˜t,ift∈BL,Y˜t=U,(5)

where [Image: pcbi.1005856.e053.jpg](W˜t)t is a Brownian motion independent of (Wt)t. Again, the mean drift across PR and BL intervals is [Image: pcbi.1005856.e054.jpg]νB*≔lb·νBlb+lp.

Fig 12View larger versionFig 12. The HBMi.The perception process P, the background process B and the resulting percept. A: A simulation on [0, 500]. B: The same realization, zoomed in on the time interval [408, 458]. Stable phases indicated by white background, unstable phases indicated by gray background. The beginnings of stable and unstable phases are marked with light blue and blue arrows, respectively.

The background process B evokes changes between the stable and the unstable state. Specifically, at the time of a percept change t*, the question of whether the process stays in the former state (S or U) or switches to the other state depends only on the value of B. Two borders, bU < bS determine this switching as follows (see Fig 12). If the former state is S, the process remains stable if and only if [Image: pcbi.1005856.e055.jpg]Bt*>b˜S (first light blue arrow in panel A), while switching to the unstable state if [Image: pcbi.1005856.e056.jpg]Bt*≤b˜S (blue arrow, panel A). Analogously, if the former state is U, the process switches to S if and only if [Image: pcbi.1005856.e057.jpg]Bt*>b˜U (right light blue arrow, panel A), while staying in U if [Image: pcbi.1005856.e058.jpg]Bt*≤b˜U (blue arrow, panel B). After the percept change, the background process B is reset to zero and then follows its usual dynamic (Eq (5)), i.e., the sign of its drift changes if and only if the state has changed. Finally, as the perception process P fluctuates between ±bS in the stable state and between ±bU in the unstable state, the value of P is reset when the state changes, to the value sgn(P)bS when changing to the stable state and to sgn(P)bU when changing to the unstable state.

Discussion of assumptions and interpretation of parameters

The technical advantage of the HBMi is that the resulting dominance times agree in most parts with the dominance times resulting from the HMM assumptions, which allows model fitting also to short data sections and comparison across clinical groups. In addition, the HBMi also provides a relation to potential underlying neuronal processes, as discussed in the following and illustrated in Fig 13.

Fig 13View larger versionFig 13. Motivation of HBMi model assumptions.For explanation and details see text.

Both HBMi-processes P and B are assumed Brownian motions with drift which may be interpreted as the activity difference between neuronal populations. Implicitly, this assumes mechanisms of self-excitation, cross-inhibition and adaptation across these neuronal populations, as proposed by various authors [14, 18]. Without explicitly modeling such mechanisms in order to reduce the number of parameters and allow model fitting, the parameter sets are reduced to the mean drifts ν and the borders b. Analogously to the HBMc, the speed of the drifts could be considered related to the inverse of the connection strengths within and across populations that engage in self excitation and cross inhibition. The border b, in analogy to the HBMc, could be considered related to the size of the respective populations under consideration. The use of different borders allows fitting of highly various response patterns and can be motivated as follows.

In the HBMi, the perception process P has two borders, bS > bU for the stable and the unstable state. This suggests different population sizes of neurons involved in the stable and unstable state. Typically bS > b > bU, suggesting that in the stable state, the activity of the dominant population is increased by joining additional neurons to the population, for example by positive feedback mediated by population S. Vice versa, in the unstable state, only a minimal population is involved in the respective percept, leading to fast changes. Thus, one could assume that the dominant percept population size is decreased by the population U (red arrows). The active population sizes are indicated by different circle sizes in the first line of Fig 13 and are assumed modulated by the background populations S and U.

The background process B models the activity difference between S and U and is also associated with two borders, [Image: pcbi.1005856.e059.jpg]b˜S and [Image: pcbi.1005856.e060.jpg]b˜U. The assumption regarding resetting of B at percept change is technically necessary to generate independent dominance times and to thus allow straightforward model fitting (cmp. section ‘Parameter estimation’). In the picture of Fig 13 it can be motivated as follows. Population S is capable of offering positive feedback to the currently active population, L or R, which results in an increased population size as described above. S is also activated by the active population. Therefore, a percept change causes a reseting to zero. However, if S had shown high previous activation (above [Image: pcbi.1005856.e061.jpg]b˜S), the activity of S can increase rapidly again, causing another stable dominance time. In contrast, in case of weak previous activation (below [Image: pcbi.1005856.e062.jpg]b˜S), the unstable population U is taking over, marking the transition to an unstable state. With opposite signs, i.e., negative drift and a small new border [Image: pcbi.1005856.e063.jpg]b˜U, the process proceeds analogously. Similar to the mean drift terms νS and νU, the drift νB is not necessarily constant but describes the mean drift of B during the period of blank display.

In addition to the potential neurophysiological interpretations of the model parameters, we give here a relation of the parameters to the response patterns. Interestingly, the seven HBMi parameters allow the reproduction of highly variable response patterns as are also observed in the empirical data sets (e.g., Fig 1). The following quantities, which are easily derived from the parameters, offer a straightforward pattern interpretation.

First, the parameter sets [Image: pcbi.1005856.e064.jpg](bS,νS*) and [Image: pcbi.1005856.e065.jpg](bU,νU*) can be interpreted analogously to the parameters (b, ν₀) in continuous stimulation. That means, an increase in the border (bS or bU) increases the mean dominance time and decreases the CV in the respective state. An increase in the drift ([Image: pcbi.1005856.e066.jpg]νS* or [Image: pcbi.1005856.e067.jpg]νU*) decreases the mean dominance time, while also decreasing the CV. Recall that the CVs of dominance times during stable and unstable states are given by

Fig 14 illustrates examples with small [Image: pcbi.1005856.e069.jpg]CVS* (panels A-D) and large [Image: pcbi.1005856.e070.jpg]CVS* (panels E-H).

Fig 14View larger versionFig 14. Impact of HBMi parameter values on the response patterns.Examples of simulated response patterns are shown for different values of the three quantities [Image: pcbi.1005856.e071]CVS*, [Image: pcbi.1005856.e072]2bSνB*/b˜SνS* and [Image: pcbi.1005856.e073]b˜U. The quantities for panels A-H were [Image: pcbi.1005856.e074]CVS*={0.2,0.2,0.2,0.2,1,1,1,1}, [Image: pcbi.1005856.e075]2bSνB*/b˜SνS*={4,4,0.4,0.4,4,4,0.4,0.4}, [Image: pcbi.1005856.e076]b˜U={0,3,0,3,0,3,0,3}.

Second, the parameters [Image: pcbi.1005856.e077.jpg]b˜S and [Image: pcbi.1005856.e078.jpg]νB* can be interpreted best when compared to bS and [Image: pcbi.1005856.e079.jpg]νS* as follows. Consider the expected fraction of [Image: pcbi.1005856.e080.jpg]b˜S reached by the background process at the end of a stable dominance time,

which is related to the transition probability from stable to unstable state. In case of a small background border [Image: pcbi.1005856.e082.jpg]b˜S<bS and small [Image: pcbi.1005856.e083.jpg]νS*, the probability of B crossing [Image: pcbi.1005856.e084.jpg]b˜S until percept change is high, such that the process remains stable. Fig 14A, 14B, 14E and 14F show such parameter combinations. An analogous term can be derived in comparison to the parameters [Image: pcbi.1005856.e085.jpg]b˜U and [Image: pcbi.1005856.e086.jpg]νU*.

Third, the parameter [Image: pcbi.1005856.e087.jpg]b˜U is related to the number of dominance times in the unstable state before changing to the stable state. Recall that the drift of B is negative during unstable periods. Therefore, a large value of [Image: pcbi.1005856.e088.jpg]b˜U implies a low probability to reach [Image: pcbi.1005856.e089.jpg]b˜U until the percept change. This implies a high expected number of dominance times in the unstable state, or a low transition probability from the unstable to the stable state. Fig 14B, 14D, 14F and 14H show examples with large values of [Image: pcbi.1005856.e090.jpg]b˜U.

Relation of the HBMi to the two state HMM

The relation of the HBMc to the one state HMM is simple as is represents only a reparameterization. Both the one state HMM and the HBMc yield independent and IG distributed dominance times. For the intermittent case, the relation of the HBMi to the two state HMM is not as straightforward. The two models are highly similar in the sense that they use two parameters to describe long and short dominance times, respectively (e.g., (μS, σS) and [Image: pcbi.1005856.e091.jpg](bS,νS*) for the stable state). In the HMM, the dominance times are IG distributed, given the state with the respective parameters. In the HBMi, the dominance times are approximately IG distributed, where the minor deviation from the IG distribution originates from the minor deviation of P from a Brownian motion with drift [Image: pcbi.1005856.e092.jpg]νS* (or [Image: pcbi.1005856.e093.jpg]νU*), instead of exactly assuming drift ν₀ during stimulation and νS (or νU) during blank display. However, the marginal distribution of P at multiples of such intervals lp + lb is identical to the marginal distribution of a Brownian motion with drift [Image: pcbi.1005856.e094.jpg]νS* (or [Image: pcbi.1005856.e095.jpg]νU*) at these time points, and the differences can only be observed in the meantime. Because dominance times usually span multiple trials of duration lp + lb, the approximation is very close. As another similarity, both models use additional parameters ((pSS, pUU) and [Image: pcbi.1005856.e096.jpg](b˜S,b˜U,νB)) to describe the transition probabilities between the stable and the unstable state.

The main difference between the HBMi and the two state HMM concerns the dynamic of the state transitions between stable and unstable state. In the HMM, transition probabilities are given by (1 − pSS) and (1 − pUU) and are independent of the duration of the previous dominance time. In contrast, in the HBMi, a transition from stable to unstable state requires that B has not reached [Image: pcbi.1005856.e097.jpg]b˜S at the end of the respective dominance time. Therefore, the transition probability [Image: pcbi.1005856.e098.jpg]p˜SU(di) depends on the duration di of the i-th dominance time, where shorter dominance times yield higher transition probabilities. Note that the position of B at the end of a dominance time di is given by an increment of a Brownian motion with drift [Image: pcbi.1005856.e099.jpg]νB* in the fixed time interval di. Therefore the position is normally distributed with mean [Image: pcbi.1005856.e100.jpg]di·νB* and variance di and the probability to remain in the stable state (which is the probability that the background process exceeds [Image: pcbi.1005856.e101.jpg]b˜S) is given by

[Image: pcbi.1005856.e102.jpg]p˜SS(di)=P(Yi+1=S|Yi=S,di)≈1-ΦνB*di,di(b˜S),(6)

where Φμ,σ²() denotes the distribution function of the normal distribution with mean μ and variance σ² and Yi is the hidden state of the i-th dominance time. Analogously, for the transition from unstable to stable state, the transition probability [Image: pcbi.1005856.e103.jpg]p˜US(di) tends to decrease with longer dominance times, and we find

[Image: pcbi.1005856.e104.jpg]p˜UU(di)=P(Yi+1=U|Yi=U,di)≈Φ-νB*di,di(b˜U).(7)

Note that we use the approximate sign ‘≈’ because the drift of B is not exactly [Image: pcbi.1005856.e105.jpg]νB* throughout, but is assumed to change between ν₀ and νB during stimulation and blank display, respectively, yielding a mean drift of [Image: pcbi.1005856.e106.jpg]νB*. Analogously to the above explanation, differences caused by the approximation can be considered minimal.

In order to obtain quantities comparable to the transition probabilities pSS and pUU in the HMM, we derive the marginal transition probabilities in the HBMi as the expected value of [Image: pcbi.1005856.e107.jpg]p˜SS and [Image: pcbi.1005856.e108.jpg]p˜UU. As shown by [39], the positions XS and XU of B at the end of an independent stable or unstable IG distributed dominance time follow the Normal Inverse Gauss (NIG) distribution. The resulting transition probabilities in the HBMi can then be calculated as

[Image: pcbi.1005856.e109.jpg]pSS*=P(XS>b˜S)andpUU*=P(XU≤b˜U).(8)

where XS is NIG-distributed with parameters (0, [Image: pcbi.1005856.e110.jpg]νS*2+νB*2, [Image: pcbi.1005856.e111.jpg]νB*, 2bS) and XU is NIG-distributed with parameters (0, [Image: pcbi.1005856.e112.jpg]νU*2+νB*2, [Image: pcbi.1005856.e113.jpg]-νB*, 2bU).

One should note that due to the difference in transition probabilities of the two models, the parameters [Image: pcbi.1005856.e114.jpg](bS,νS*) are not direct reparameterizations of (μS, σS) (and similarly for the unstable state). Furthermore, the dependence of the transition probability on the length of the previous dominance time is one important new aspect of the HBMi not described in the HMM, which will also be used in the section ‘Application of the HBM to the sample data set’ for comparison of models and empirical observations.

Parameter estimation

In order to estimate the HBMi parameter set [Image: pcbi.1005856.e115.jpg]Θ=(bS,νS*,bU,νU*,b˜S,b˜U,νB*), we use ML estimation. The likelihood L is given by

[Image: pcbi.1005856.e116.jpg]L(d₁, …, dn|Θ) ≈ αS(n) + αU(n), (9)

where the forward variable αj(i) denotes the probability of being in state j at time i while observing (d₁, …, di). The approximation is again due to the averaging of drifts during blank display and stimulus presentation. The forward variables can be derived recursively [34, 40] by

[Image: pcbi.1005856.e117.jpg]αj(i)={πjfj(di)ifi=1fj(di)∑j∈{S,U}αj(i-1)p˜ij(di-1)ifi>1(10)

with πj as starting distribution and fS and fU denoting the densities of [Image: pcbi.1005856.e118.jpg]IG(2bS/νS*,2bS/νS*3) and [Image: pcbi.1005856.e119.jpg]IG(2bU/νU*,2bU/νU*3) distributed random variables, respectively. Details on the maximization algorithm can be found in the Materials and methods. After estimation of Θ, the estimates of [Image: pcbi.1005856.e120.jpg]ν^S and [Image: pcbi.1005856.e121.jpg]ν^U can be obtained using Eq (4), and analogously for [Image: pcbi.1005856.e122.jpg]ν^B.

Precision of parameter estimation. The variability of the parameter estimates in the HBMi is studied analogously to the HMM, using the RE and the AE of the parameter estimates obtained in 1000 simulations of the 61 parameter combinations estimated from the empirical data set. For the border parameters bS, bU, [Image: pcbi.1005856.e123.jpg]b˜S, [Image: pcbi.1005856.e124.jpg]b˜U, we use the RE, while for the typically small drift parameters [Image: pcbi.1005856.e125.jpg]νS*, [Image: pcbi.1005856.e126.jpg]νU*, [Image: pcbi.1005856.e127.jpg]νB*, the AE is used. Again, one set of simulations uses the time horizon of the empirical data, T₁ = 1200s (Fig 15A), and a second simulation was performed using T₂ = 3600s (Fig 15B). According to the simulation results, the precision of parameter estimation in the given parameter range is not always satisfactory for the given parameterization [Image: pcbi.1005856.e128.jpg](bS,bU,b˜S,b˜U,νS*,νU*,νB*), yielding average errors (i.e., mean REs or AEs across all variables) smaller than 0.25 in only 24 out of 61 cases for the empirical time horizon T₁ and still in only 44 out of 61 for the tripled sample size of T₂. This suggests that these raw parameters yield less reliable estimates because different combinations of b and ν can yield the same mean stable dominance time. In contrast, the following set of derived parameters that are more easily comparable to the HMM parameters shows better properties. We denote by [Image: pcbi.1005856.e129.jpg](μS*,σS*) and [Image: pcbi.1005856.e130.jpg](μU*,σU*) the mean and standard deviation of dominance times in the HBMi in the stable and unstable state, respectively. These have analogous interpretations as (μS, σS) and (μU, σU) in the HMM and are given by

[Image: pcbi.1005856.e131.jpg]μS*≈2bS/νS*andσS*≈2bS/νS*3,(11)

and analogously for [Image: pcbi.1005856.e132.jpg]μU* and [Image: pcbi.1005856.e133.jpg]σU*, where the approximation is again due to the minimal difference between the mean drift [Image: pcbi.1005856.e134.jpg]νS* and the changing drift νS + ν₀ during presentation and blank display. In addition, we consider the transition probabilities between the states, [Image: pcbi.1005856.e135.jpg]pSS* and [Image: pcbi.1005856.e136.jpg]pUU* (cmp. Eq (8)). Fig 15C and 15D show the median REs for [Image: pcbi.1005856.e137.jpg]μS*, [Image: pcbi.1005856.e138.jpg]σS*, [Image: pcbi.1005856.e139.jpg]μU*, [Image: pcbi.1005856.e140.jpg]σU* and median AEs for [Image: pcbi.1005856.e141.jpg]pSS*, [Image: pcbi.1005856.e142.jpg]pUU* obtained in the 1000 simulations of length T₁ (panel C) and T₂ (panel D). Concerning this parameterization, 51 parameter combinations yielded average errors smaller than 0.25 across all variables for T₁, while as many as 58 out of 61 cases showed errors smaller than 0.25 for T₂. These simulations suggest that the parameterization [Image: pcbi.1005856.e143.jpg](μS*,σS*,μU*,σU*,pSS*,pUU*) yields more reliable parameter estimates in the HBMi than the original parameterization of borders and drifts.

Fig 15View larger versionFig 15. Precision of parameter estimates in the HBMi.For each of the 61 parameter constellations estimated from the sample data set, 1000 simulations were performed with the HBMi. (A) and (B): log(median(REs)) of the border parameters bS, bU, [Image: pcbi.1005856.e144]b˜S, [Image: pcbi.1005856.e145]b˜U, and log(median(AEs)) of the drift parameters [Image: pcbi.1005856.e146]νS*, [Image: pcbi.1005856.e147]νU*, [Image: pcbi.1005856.e148]νB* with T = 1200s (A) and T = 3600s (B). (C) and (D) log(median(REs)) for the derived model parameters [Image: pcbi.1005856.e149]μS*, [Image: pcbi.1005856.e150]σS*, [Image: pcbi.1005856.e151]μU*, [Image: pcbi.1005856.e152]σU* and log(median(AEs)) for the parameters [Image: pcbi.1005856.e153]pSS* and [Image: pcbi.1005856.e154]pUU* for T = 1200s (C) and T = 3600s (D). Parameter combinations with mean errors < 0.25 across all variables plotted in black.

Application of the HBM to the sample data set

For the analysis of the sample data set, the HBMc and HBMi were fitted to the dataset of [10] containing responses to continuous and intermittent stimulation for each of 29 patients with schizophrenia and 32 control subjects, using the parameter estimation described in the section ‘Parameter estimation’.

Model fit. Because the HBMc represents only a reparameterization of the one-state HMM, results are completely analogous for continuous presentation. Thus, a high variability of response patterns can be described with the two parametric distribution (see Fig 16A and 16B), including also different means and variances of dominance times. For intermittent presentation, the HBMi and the two-state HMM are similar, but also show a number of differences. As a first similarity to the two-state HMM, the HBMi can also describe and reproduce a high variety of response pattens (Fig 16C–16F). For example, these include highly regular stable states that may or may not be interrupted by short unstable phases (C,E) or response patterns with different degrees of regularity (D, F) and different alternation rates. Note also that the response patterns to intermittent stimulation of six out of the 61 subjects were described better by the one-parametric HBMc (e.g., E).

Fig 16View larger versionFig 16. HBM: Examples of simulated response patterns.Example empirical response patterns from Fig 1 and response patterns simulated by the HBM (orange). Responses to continuous and intermittent presentation are plotted in green and blue, respectively. The estimated parameters used for simulation are given in Tables 3 and 4, respectively. (A) and (B) Continuous presentation and HBMc. (C)-(F): Intermittent presentation and HBMi. Below the orange response patterns, one can see the perception process P and (for the HBMi) the background process B corresponding to the respective simulation.

See full table Table 3. Estimated HBMi parameter combinations of the typical response pattens to continuous presentation shown in Fig 1A and 1B.

See full table Table 4. Estimated HBMi parameter combinations of the typical response pattens to intermittent presentation shown in Fig 1C–1F.

In addition to the close description and reproduction of the patterns in the empirical data, one interesting additional aspect is captured by the HBMi, which cannot be described in the two-state HMM. As explained in the section ‘Relation of the HBMi to the two state HMM’, the probability of a transition from stable to unstable state decreases with the length of the dominance time in the HBMi. The same is true for the reverse transition. Indeed, the same observation can be made in the empirical data set, while this dependence cannot be captured within the HMM (Fig 17).

Fig 17View larger versionFig 17. Mean dominance times in stable state as a function of the successive state.In the empirical data set, the mean dominance time before a state transition SU is shorter than the mean dominance time before SS (A). This observation can be reproduced in the HBMi (C), while the expected dominance time in the HMM (B) is equally long for intervals with and without transition, i.e., independent from transition between states. In the empirical data set [10], the state transitions were estimated using the Viterbi paths [34]. Analogous results were obtained using a fixed threshold. For the HBMi, the expected dominance times were derived from the HBMi parameters according to Eq (15) in the Materials and methods. In the HMM, expected dominance times correspond to [Image: pcbi.1005856.e165]μ^S.

Group differences. Due to the high correspondence of the HBMi response patterns with the empirical data and the neurophysiologically related parameterization, the HBMi provides potential additional links to underlying neuronal processes of the observed differences between control subjects and patients with schizophrenia. Here, we consider three aspects related to continuous and intermittent stimulation and to the transition between these two conditions.

First, concerning continuous stimulation, we note that the one-state-HMM and the HBMc yield identically distributed sequences of dominance times. Note that the mean dominance time [Image: pcbi.1005856.e166.jpg]μ^ in the HMM (Fig 9A) therefore equals the corresponding value [Image: pcbi.1005856.e167.jpg]2b^/ν^0 in the HBMc. As a consequence, the results and interpretation were identical, i.e., the higher alternation rate of the control subjects during continuous presentation (Fig 2) was reflected in a smaller value of [Image: pcbi.1005856.e168.jpg]2b^/ν^0. When analyzing the individual parameters b and ν₀, we found no group differences in the drift ν₀, but a tendency for a larger neuronal pool b involved in sensory processing in the group of patients with schizophrenia (Fig 18A, p < .1, two-sided Wilcoxon test).

Fig 18View larger versionFig 18. Group differences between patients with schizophrenia and control subjects in the HBM.The HBMc parameter [Image: pcbi.1005856.e169]b^ (A) and the relation [Image: pcbi.1005856.e170]b^/b^S (B). (C): Differences in [Image: pcbi.1005856.e171]b^S/b˜^S and [Image: pcbi.1005856.e172]p^SS*. The raw data together with the median (colored diamonds) and 25%- and 75%-quantiles are shown. All p-values of a two-sided Wilcoxon test were below.1.

Second, regarding intermittent presentation, we focussed on the derived parameters [Image: pcbi.1005856.e173.jpg](μS*,σS*,μU*,σU*,pSS*,pUU*), which show correspondence to the HMM and good estimation properties, instead of testing the border and drift parameters individually. Similar to the HMM, we found an increased mean unstable dominance time in the group of patients with schizophrenia as compared to healthy controls. More importantly and also consistent with the HMM, we also found that patients with schizophrenia showed a decreased relative time spent in S, [Image: pcbi.1005856.e174.jpg]φS* (p < .1, two-sided Wilcoxon test, see Eq (17) in the Materials and methods). Again, the probability [Image: pcbi.1005856.e175.jpg]pSS* to stay in the stable state seemed to be the main variable contributing to this difference, being significantly reduced in patients with schizophrenia (Fig 18C, p < .1, two-sided Wilcoxon test). In order to identify potential neurophysiological mechanisms underlying this group difference, we further investigated the components of [Image: pcbi.1005856.e176.jpg]pSS*. Noting that no difference was observed in the mean stable dominance time [Image: pcbi.1005856.e177.jpg]μS*=2bS/νS* and in the relation [Image: pcbi.1005856.e178.jpg]νS*/νB* of the drift in stable state and the drift of the background process, one interesting parameter is [Image: pcbi.1005856.e179.jpg]bS/b˜S. Keeping all other parameters constant, an increase in this quantity means an increase in [Image: pcbi.1005856.e180.jpg]pSS* because the crossing of [Image: pcbi.1005856.e181.jpg]b˜S at the end of a stable dominance time, i.e., staying in the stable state, gets more likely. In the empirical data set, we found increased values of [Image: pcbi.1005856.e182.jpg]bS/b˜S in the control group (Fig 18C, p < .1, two-sided Wilcoxon test). In terms of the potential neurophysiological interpretation, this would suggest an increased population size involved in stable perception processing as a potential main underlying mechanism.

Third, a similar observation resulted from the derived parameter bS/b, which describes a transition of population sizes from continuous to intermittent stimulus processing. In the data set, we observed increased values of bS/b for the control group (Fig 18B) as compared to the group of patients with schizophrenia. Again, applying a potential neurophysiological interpretation, this would suggest that in the control group, the neuronal pool that is additionally recruited during intermittent presentation in stable states could be higher than in the group of patients with schizophrenia. Note that this result could not be obtained in the HMM, which does not explicitly describe a potential transition between continuous and intermittent presentation. However, note that this result should be interpreted carefully due to reduced estimation precision of the parameters [Image: pcbi.1005856.e183.jpg]bS/b˜S and bS/b.

In summary, this analysis of the HBM model parameters suggests the following explanation for the observed phenomenon that the alternation rate of the perceived percepts is increased for patients with schizophrenia during intermittent stimulation, while being decreased during continuous stimulation. In general, the relative time spent in the unstable state was increased in patients with schizophrenia. According to the HBM, this was attributed to an increased probability of transition from the stable to the unstable state, which could be potentially related to a decreased recruitment of neurons in the stable state. More specifically, a larger neuronal pool is hypothesized to account for the increased stability, and we accordingly observe a smaller increase in the neuronal pool from continuous to intermittent stable presentation in the patients with schizophrenia, which is suggested by the HBM as a main mechanism for the observed group differences. This analysis of a potential underlying mechanism, explaining the observed group differences also in the transition from continuous to intermittent presentation, is a particular advantage of the HBM over the HMM, because the HBM provides mechanistic explanations and variables with potential neurophysiological interpretations.

Discussion

Summary and implications

In the present article we have proposed a model framework for the description and analysis of perceptual responses to bistable stimuli. In particular, the first goal was to describe a high number of observed patterns in responses to continuous and intermittent stimulation and their differences between a group of patients with schizophrenia and healthy controls. The variety of patterns includes more or less regular dominance times during continuous stimulation and a switching between long and short dominance times, i.e., stable and unstable states, during intermittent stimulation, with a tendency for periodically occurring percept changes.

We started on a descriptive level, assuming that dominance times were generated by a simple HMM with only one state for continuous presentation and a stable and an unstable state in intermittent presentation. The HMM was sufficiently small to allow model fit to short empirical data sets and could also describe the high variety of empirically observed response patterns in continuous and intermittent presentation. Interestingly, it also revealed a high degree of reproducibility of response patterns of the same subject across different sessions. In addition, it allowed to relate observed group differences in the rate of percept alternations to HMM parameters, suggesting that especially the relative time spent in the stable state was reduced in the patients with schizophrenia.

Our second goal was to relate the observed response patterns and group differences to potential underlying mechanisms and thus, to build a link to models with detailed neurophysiological assumptions [7, 13, 14, 18, 20] that may not include all types or response patterns and/or may not allow fitting to short data sets. To that end, we proposed a hierarchical model of interacting Brownian motions (HBM). The HBM is based on the common assumption that the sequence of percept changes results from a competition of conflicting neuronal populations [14, 17, 18, 20, 27]. Instead of modeling these in detail, we describe the activity difference by a Brownian motion P with drift ν₀ [21] between two borders ±b, where the first hitting times of the borders indicate percept changes. Roughly speaking, the drift ν₀ could be considered related to the neuronal interactions within and between the populations, while the border b could be considered related to the population sizes. In order to describe responses to intermittent presentation, this mechanism is adapted in another population pair. These populations exhibit a corresponding background process B that evokes switching between stable and unstable states, similar to switching between the two percepts. In particular, B causes the perception process P to change parameters from small drift νS and large border bS in the stable state to fast drift νU and small border bU in the unstable state.

The HBM could be fitted nicely to the given empirical data set, reproducing a high variety of response patterns to continuous and intermittent stimulation in healthy subjects and patients with schizophrenia. In particular, the model fit was even improved over the descriptive HMM by reproducing shorter stable dominance times before a change to the unstable state. The HBM also provided more detailed explanations for the observed group difference that patients with schizophrenia showed higher alternation rates during intermittent stimulation, while percept alternation was decreases during continuous presentation. In particular, the HBM contains additional mechanisms of switching between stable and unstable state for intermittent presentation, which is assumed inactive during continuous presentation. The HBM, similar to the HMM, suggests an increased probability of switching to the unstable state for the patients with schizophrenia and thus, a longer relative time spent in the unstable state. The HBM also provides additional potential explanations related to the borders, or assumed population sizes, suggesting a higher increase from continuous (border b) to stable intermittent presentation (border bS) in the healthy subjects. This is a first finding on the transition from continuous to intermittent presentation, which results from including both continuous and intermittent presentation in one model.

These findings suggested by the HBM, which include a longer relative time spent in the unstable state for the patients with schizophrenia and a smaller population size involved in percept stabilization, are also in agreement with recent findings of [41]. They studied the learning behavior of healthy subjects of whom the degree of delusional ideation [42] had been measured. In compliance with earlier studies [for a review see 43], they reported that subjects with larger delusion proneness made decisions on the basis of less information and were also less resilient against irrelevant information [compare also the literature about jumping to conclusions, e.g. 44, 45]. In the present setting, the ambiguous stimulus represents a constant source of partly contradicting visual information [see also 5, 8]. In that sense, the unstable state could be considered a state in which one is less resilient against this contradicting visual information, which yields a high rate of percept changes. The fact that the patients with schizophrenia spent more time in the unstable state is therefore highly consistent with the findings of [41]. Moreover, this finding is also compatible with current models of schizophrenia in the framework of predictive coding [46] that propose a reduced top-down influence of stored predictions. However, it goes beyond previous work by highlighting the role of a background process that controls the balance between stable and unstable states in perceptual inference. In addition, the population sizes could be considered related to the amount of information taken into consideration to create a percept. Again, consistently with [41], we find, in the stable state, larger estimated population sizes, bS, of the perceptual populations L and R in healthy controls than in patients with schizophrenia. Also, these population sizes are typically much larger than the population sizes in the unstable state (bS >> bU), which would be consistent with the notion that subjects in the unstable state need less information to change their perception.

Applicability and model extensions

The HBM may also be used to describe dominance times resulting from other experiments with ambiguous visual stimuli studying, e.g., motion-induced-blindness, binocular rivalry, moving plaids, the Necker Cube, orthogonal gratings or the house/face-paradoxon [e.g. 21, 47] or also bistable auditory stimuli [48]. The HBM is, however, not designed for tristable stimuli, and transient stimulus manipulations as used in after-effect studies cannot be captured by the HBM in its current form. In different bistable settings, the HBM cannot be applied directly, but would allow for potential extensions. For example, in its present form, the HBM describes only balanced perception. However, it could be extended with respect to unbalanced bistable displays, e.g., for different eye contrasts during binocular rivalry [49], by choosing different drift parameters for the positive and the negative drift direction during presentation. Similarly, the drift could be chosen to vary as a function of attention [50, 51] or as a function of long-term history (e.g., the cumulative history H proposed in [31]). In studies on mixed perception during binocular rivalry [19], one might use an additional border to define an intermediate range for the perception process in which mixed perception is described.

One should note that the HBMi in its current form is restricted to a duration of blank displays lb ≤ lp ⋅ ν₀/νS. For longer blank displays, the mean drift of P during stable states, [Image: pcbi.1005856.e184.jpg]νS*, will be negative, yielding no perception change with high probability. However, it would be possible to extend the model accordingly, assuming a temporal evolution in the drift parameters, given corresponding extended empirical observations. In addition, note that the border of the perception process is assumed to be b during continuous stimulation and bS (or bU) during intermittent stimulation. Therefore, an instantaneous change of intermittent to continuous presentation is not yet described. Here, we qualitatively assume that the border jumps very fast from b to bS with the onset of a blank display, while going back slowly during stimulation. A transition from continuous to intermittent presentation would therefore instantly change the response pattern, while a reverse transition would gradually reverse the change back to the one-state process. Quantitative validation and fitting of this assumption would be interesting, but requires corresponding empirical observations, in which the length of the presentation period lp is varied. This would also allow investigation of potential relations between the HBMc and HBMi parameters and thus, between the mechanisms assumed to underlie the identified group differences.

Concerning the impact of the duration of the blank display lb, two aspects should be discussed. First, the HBM can theoretically reproduce a phenomenon reported earlier in [14]. Conditional that one percept has been present for a short while, the probability of a percept change rises with the blank duration lb. In the HBMi, the same is observed during the unstable state with typically short dominance times: During the unstable state the drift in the blank displays, νU, is typically larger than the drift ν₀ during stimulation. Therefore, longer blank displays speed up P, thereby reducing perceptual stability.

Second, one interesting potential model extension is concerned with the relationship between the length of the blank display and the alternation rate. As reported earlier by [14, 15, 18], the mean dominance time in intermittent presentation has been found to be a function of the relationship between the presentation length lp (or ‘ON’-period) and the length of the blank display lb (or ‘OFF’-period). Particularly, the dependence between lb and the alternation rate is non-monotonic, as would be implied in the HBMi, but follows an inverted U-shape [22, 24, 25] with a peak roughly at 0.4 s. Such an inverted U-shape would be possible in a model extension of the HBMi. As discussed in the Results, the drift terms νS, νU only represent the mean drift across the period of blank display, which is sufficient and parsimonious in the given data set with fixed length of blank display. However, the model would be fully consistent with the assumption that the drifts change during the ‘OFF’-period, such that the mean drifts νS(lb) and νU(lb) are functions of the length of the blank display lb. In Fig 19A these mean drifts νS, νU decrease with lb, where the stronger drifts at the beginning of the blank display could be effects of the recent stimulation. Panel B shows the resulting mean alternation rate, which has an inverted U-shape with a maximum around 0.4 s and shows increased stability under intermittent stimulation for lb > 0.7. This increased stability is caused first by a small drift νS < ν₀ in that range. Second, it is also caused by the fact that the time interval lb in which the background process B has positive drift is longer, leading to an increased probability to reach [Image: pcbi.1005856.e185.jpg]b˜S and thus, to stay in the stable state. Estimation of the functions νS(lb) and νU(lb) from a suitable data set with variable lengths of blank displays would be an interesting task.

Fig 19View larger versionFig 19. HBMi extension for different lengths of blank display.Parameters derived for subject B and C in Tables 3 and 4, which is the same subject during continuous and intermittent presentation. A) The mean drifts νS, νU decrease with the length of the blank display lb. B) The resulting mean alternation rate ρ under intermittent stimulation (black) is derived using formula (18) (Materials and methods). The empirical mean alternation rates per minute of the subject during intermittent viewing with lb = 0.8 and during continuous viewing are marked by blue and green lines, respectively.

In summary, the proposed HBM intends to provide a link between empirical data analysis and mechanistic modeling. On the one hand, it aims at precisely describing the high variety of response patterns observed in perceptual responses to bistable stimuli. On the other hand, it aims at bridging the gap to detailed mechanistic models of bistable perception, allowing assumed processes to be fitted to short empirical data sets and thus, also the analysis of group differences. Various extension possibilities show a potential of the HBM to investigate related experimental contexts. By including both continuous and intermittent stimulation, the HBM can thus provide interesting new hypotheses on potential neuronal mechanisms of cognitive phenomena.

Materials and methods

Estimation of HMM parameters

Here we describe the estimation procedures of the HMM parameters for continuous presentation and for intermittent presentation. We denote by d ≔ (d₁, d₂, …, dn) the set of dominance times modeled as realizations of random variables D = (D₁, …, Dn).

Continuous presentation: ML estimation

The HMM for continuous presentation assumes that all dominance times di are independent and Inverse Gaussian distributed with mean μ and standard deviation σ. The log-likelihood function is then given by

[Image: pcbi.1005856.e186.jpg]logL(d|μ,σ)=12∑i=1nlog(μ32πσ2di3)-μσ2(di-μ)2di.(12)

Setting the partial derivatives with respect to μ and σ to zero yields the estimates of Eq (12) [32]:

Intermittent presentation: Baum-Welch-Algorithm

The HMM for intermittent presentation uses the parameter set ϑ = (μS, σS, μU, σU, pSS, pUU) for the mean and standard deviation in the stable and unstable state, respectively, and the transition probabilities between these two states. This parameter set is estimated with the Baum-Welch-Algorithm [BWA, 33] which is an iteratively working instance of the EM-Algorithm maximizing the model likelihood locally. Here, we explain it briefly [for details see, e.g., 34]. In the first step one applies the so called Forward- and Backward-Algorithm. The forward-variable αj(i) is defined as the probability of observing the sequence d₁, d₂, …, di and being in state j at time i, given the model parameters. The backward-variable βj(i) denotes the analogous probability of observing the ending partial sequence di+1, di+2, …, dn and being in state j at time i. To avoid underflow we normalize both variables [e.g. 34], resulting in the normalized variables [Image: pcbi.1005856.e188.jpg]α˜j(i), [Image: pcbi.1005856.e189.jpg]β˜j(i). The normalized variables can be derived iteratively as follows

where πj denotes the probability that the Markov chain starts in state j, pSU = 1 − pSS, pUS = 1 − pUU and [Image: pcbi.1005856.e192.jpg]fμ,σIG(x) denotes the density of the IG distribution with expectation μ and standard deviation σ evaluated at x. Note that we suppress the dependence of [Image: pcbi.1005856.e193.jpg]α˜j(i) and [Image: pcbi.1005856.e194.jpg]β˜j(i) on the parameter set ϑ for convenience.

The forward and backward variables are used to derive the probability γj(i|ϑ) of being in state j at time i, given the whole sequence d ≔ (d₁, …, dn) and the parameters ϑ

Moreover, we need the probability ξj,k(i|ϑ) of being in state j at time i and in state k at time i + 1, given the data d and the parameters ϑ,

To iteratively derive the parameter estimates, the BWA applies expectation maximization as follows. Let ϑ(m) denote the parameter estimates after the m-th iteration step, and let [Image: pcbi.1005856.e197.jpg]𝔜 denote the set of all possible state sequences of the hidden Markov chain. Let Y = (Y₁, …, Yn) denote a [Image: pcbi.1005856.e198.jpg]𝔜-valued random variable and y = (y₁, …, yn) a realization of Y. The BWA then iteratively maximizes the Q-function [e.g. 35] over Y,

i.e., the expectation of the complete-data log-likelihood L across all possible paths [Image: pcbi.1005856.e200.jpg]y ∈ 𝔜. The updated parameter set ϑ(m+1) is chosen such that it maximizes Q. For a fixed state sequence y = (y₁, …, yn) the log-likelihood of the data is

Insertion into Q yields

Note that the first line depends only on the initial distribution π, the second line depends on the transition probabilities and the third line depends on the parameters of the Inverse Gaussian distributions. Therefore, iterative parameter estimation separately maximizes these terms. Note further that we can rewrite [Image: pcbi.1005856.e203.jpg]ℙ(Yi = yi|d, ϑ(m)) = γyi(i|ϑm) and [Image: pcbi.1005856.e204.jpg]ℙ(Yi-1 = yi-1, Yi = yi|d, ϑ(m)) = ξyi-1yi(i - 1|ϑm), which yields the following estimates in the m + 1-st iteration step.

For the entries of the transition matrix, Lagrange maximization of the second line under the constraints pSS + pSU = pUU + pUS = 1 yields the estimate

in the (m + 1)-th step of the BWA [e.g. 34].

The big bracket in the last line of the Q-function is analogous to the likelihood function of the IG distribution (Eq (12)) given in the section ‘Continuous presentation: ML estimation’ (with the additional indices yi and the weighting factors [Image: pcbi.1005856.e206.jpg]ℙ(Yi = yi|d, ϑ(m))). Therefore, setting the partial derivatives to zero yields for j ∈ {S, U} the analogous estimates

which are the updates for parameters of the Inverse Gaussian distributions in the (m + 1)-th step of the BWA.

In order to update the starting distribution, we do not maximize the respective summand of Q but we assume the stationary distribution as the starting distribution π. The stationary distribution is derived using the eigenvectors of the transpose of the transition matrix [52], which yields

These iterative steps are repeated until a desired level of convergence is reached.

Starting values and constraints. As starting values [Image: pcbi.1005856.e209.jpg]μS(s), [Image: pcbi.1005856.e210.jpg]σS(s), [Image: pcbi.1005856.e211.jpg]μU(s), [Image: pcbi.1005856.e212.jpg]σU(s), [Image: pcbi.1005856.e213.jpg]pSS(s), [Image: pcbi.1005856.e214.jpg]pUU(s) for the Baum-Welch algorithm we chose, in correspondence with the data set, [Image: pcbi.1005856.e215.jpg]pSS(s)=pUU(s)=0.5; [Image: pcbi.1005856.e216.jpg]μS(s)=4; [Image: pcbi.1005856.e217.jpg]σU(s)=5. In order to reduce the probability that the Baum-Welch-Algorithm will be captured in a local extremum, we chose ten equidistant values for [Image: pcbi.1005856.e218.jpg]μS(s) ranging between 60 and 0.95maxi di, and for each value of [Image: pcbi.1005856.e219.jpg]μS(s) we choose ten equidistant values for [Image: pcbi.1005856.e220.jpg]σS(s) between 10 and [Image: pcbi.1005856.e221.jpg]1.1μS(s). Out of the resulting one hundred sets of parameter estimates we chose the parameter set with the highest log-likelihood. If the response pattern shows only dominance times larger than 30 seconds we reduce the model to the stable phase. The parameters μS and σS are derived by ML as described in the section ‘Continuous presentation: ML estimation’, and we set pSS ≔ 1. If the dominance time are only smaller than 30 seconds, we only estimate μU and σU by ML and use pUU = 1.

For subjects with relatively clear distinction between long and short dominance times this procedure yields reasonable estimates. For subjects with less clear distinction, we added the following constraints based on the idea that short dominance times should not affect estimation of the stable parameters and long dominance times should not affect estimation of unstable parameters. Note that in continuous presentation where no state exists, about 90% of the dominance times are shorter than 15 seconds, while only about two percent are larger than 30 seconds. Therefore, we first require [Image: pcbi.1005856.e222.jpg]σ^S>1 such that not just the largest dominance time is estimated as stable and all others are categorized as unstable (which may increase the likelihood). Second, we do not accept HMMs with [Image: pcbi.1005856.e223.jpg]μ^S<0.98μ^15 where [Image: pcbi.1005856.e224.jpg]μ^15≔(1/|d1d>15|)∑i=1ndi1di>15. This prevents dominance times smaller than 15 seconds to be considered for the estimation of μS. Third, we require [Image: pcbi.1005856.e225.jpg]μ^S<1.02μ^75 where [Image: pcbi.1005856.e226.jpg]μ^75≔(1/|d1d>75|)∑i=1ndi1di>75 if any dominance time is larger than 75 seconds and [Image: pcbi.1005856.e227.jpg]μ^75=75 otherwise such that rather stable dominance times longer than 75 seconds are not classified as unstable.

Expected relative time spent in the stable state

Here we derive the formula for the expected time φS spent in the stable state investigated in Fig 9. To that end, let Nj, j ∈ {S, U}, denote the (random) number of dominance times in state j before a state change. As the number of dominance times before a state change is geometrically distributed with probability 1 − pjj, we have for the expectation

Moreover, let [Image: pcbi.1005856.e229.jpg](DiS)i≥1 be a sequence of independent IG(μS, σS)-distributed random variables and [Image: pcbi.1005856.e230.jpg](DiU)i≥1 be a sequence of IG(μU, σU)-distributed random variables. We derive φS as

The length of a stable phase is a random variable distributed like [Image: pcbi.1005856.e232.jpg]∑i=1NSDiS, where [Image: pcbi.1005856.e233.jpg]D1S,…,DNSS are independent from each other and also independent of NS. Therefore we have [53]

and analogously for the expected length of an unstable phase. This yields

[Image: pcbi.1005856.e235.jpg]φS=μS/(1-pSS)μS/(1-pSS)+μU/(1-pUU)=(1-pUU)μS(1-pUU)μS+(1-pSS)μU.(13)

Estimation of HBM parameters

Here we describe the estimation procedures of the HBM parameters for continuous presentation and for intermittent presentation.

HBMc

Recall that the HBMc describes the perception process P as a Brownian motion with drift ±ν₀, where the sign of the drift changes at the first hitting time of the borders ±b. The two parameters (b, ν₀) are estimated using the ML method. Due to transformation invariance of the ML estimates, we can therefore use the ML estimators from the section ‘Continuous presentation: ML estimation’ applying that the resulting dominance times are IG distributed with expectation μ = 2b/ν₀ and variance [Image: pcbi.1005856.e236.jpg]σ2=2b/ν03, or conversely, [Image: pcbi.1005856.e237.jpg]b=(1/2)·μ3/σ2 and [Image: pcbi.1005856.e238.jpg]ν^0=μ/σ2. This yields

HBMi

The parameters of the HBMi are [Image: pcbi.1005856.e240.jpg]Θ=(bS,νS*,bU,νU*,b˜S,b˜U,νB*) where [Image: pcbi.1005856.e241.jpg](bS,νS*) and [Image: pcbi.1005856.e242.jpg](bU,νU*) denote the border and drift parameters of the perception process P in the stable and the unstable state, respectively. [Image: pcbi.1005856.e243.jpg]b˜S, [Image: pcbi.1005856.e244.jpg]b˜U, [Image: pcbi.1005856.e245.jpg]νB* are the border and drift parameters of the background process B which determines the hidden state. Recall that the likelihood of the whole model given the data and the parameter vector Θ is approximately

[Image: pcbi.1005856.e246.jpg]L(d|Θ) ≈ αS(n) + αU(n), 

using the forward variables αj(i) defined recursively in Eq (10). Note again that we suppress dependence on the parameters Θ for convenience. In practice we need to avoid underflow when calculating αj(i). To that end the forward variables are normalized such that ∑j αj(i) = 1 for all time points i, using the following steps [40]. For j ∈ {S, U},

The likelihood then derives as

[Image: pcbi.1005856.e248.jpg]L(d|Θ)≈∏i=1nci(α˜S(n)+α˜U(n))=∏i=1nci,yieldinglogL≈∑i=1nlog(ci).(14)

Parameter estimation is then obtained by maximizing ∑i log(ci), which is a function of the model parameters Θ. To that end we apply the Newton-type algorithm [54] implemented in the R-function nlm(). Alternatively the COBYLA-Algorithm [55] for maximization under non-linear constraints can be applied.

Next we discuss the set of starting values [Image: pcbi.1005856.e249.jpg]{bS(s),νS*(s),bU(s),νU*(s),b˜S(s),b˜U(s),νB*(s)} for the optimization algorithm. Let [Image: pcbi.1005856.e250.jpg]U≔{μ^15,μ^30} and [Image: pcbi.1005856.e251.jpg]O≔{σ^15,σ^30,1.15σ^30}, where [Image: pcbi.1005856.e252.jpg]μ^k and [Image: pcbi.1005856.e253.jpg]σ^k denote the empirical mean and standard deviation of all dominance times larger than k seconds. Then we choose the initial values for bS and [Image: pcbi.1005856.e254.jpg]νS* from the sets

Depending on [Image: pcbi.1005856.e256.jpg]νS*(s) and [Image: pcbi.1005856.e257.jpg]bS(s), we choose [Image: pcbi.1005856.e258.jpg]bU(s)∈{0.01bS(s),0.05bS(s),0.15bS(s)}, [Image: pcbi.1005856.e259.jpg]b˜S(s)∈{bS(s),10bS(s)}, [Image: pcbi.1005856.e260.jpg]νU*(s)∈{1.01νS*(s),3νS*(s),7νS*(s)}, furthermore [Image: pcbi.1005856.e261.jpg]νB*(s)∈{0.1,3} and [Image: pcbi.1005856.e262.jpg]b˜U(s)∈{-3,0,3}.

The initial value for [Image: pcbi.1005856.e263.jpg]πS = ℙ(process starts in the stable state) is set to [Image: pcbi.1005856.e264.jpg]π^S=1 if d₁ ≥ 45, [Image: pcbi.1005856.e265.jpg]π^S=0 if d₁ ≤ 15, and [Image: pcbi.1005856.e266.jpg]π^S=1/2 otherwise. Alternatively and comparable to the HMM the stationary distribution can be used for the initial distribution, thereby reducing the number of parameters by one. As this, however, requires numerical integration and increases the computational effort considerably and the differences between the two approaches were negligible we used the simpler rule.

The maximization algorithm is applied using all combinations of starting values. We then take the set of parameter estimates which yields the highest log-likelihood and fulfills the following constraints

Constraints A) to D) result from the model assumptions, E) prevents dominance times smaller than 15 seconds to be considered for the estimation of μS, analogously to the HMM procedure. Constraint F) prevents too big estimates of standard deviations, which would yield implausible results.

In cases in which only dominance times larger than 30 seconds are observed, we apply the algorithm described for the HBMc, where bS, [Image: pcbi.1005856.e268.jpg]νS* are estimated like in the section ‘HBMc’ and [Image: pcbi.1005856.e269.jpg]b˜S is set to zero. In cases in which only dominance times up to 30 seconds are observed, we proceed analogously, where bU, [Image: pcbi.1005856.e270.jpg]νU* are estimated like in the section ‘HBMc’ and set [Image: pcbi.1005856.e271.jpg]b˜U=1010. In either case we set [Image: pcbi.1005856.e272.jpg]νb*=10 and do not estimate the other variables.

HBMi: Dominance time before state changes

Stable dominance time before a change to the unstable state

In the section ‘Relation of the HBMi to the two state HMM’ we discussed that in the HBMi, stable dominance times are shorter when occurring directly before a state change to the unstable state than when followed by another stable dominance time. The same observation was also made in the sample dataset (Fig 17). Therefore, we derive here the expected length [Image: pcbi.1005856.e273.jpg]μS- of a stable dominance time before a state change and its corresponding expected length [Image: pcbi.1005856.e274.jpg]μS+ before another stable dominance time.

Let [Image: pcbi.1005856.e275.jpg]Θ=(bS,νS*,bU,νU*,b˜S,b˜U,νB*) be the parameter set of an HBMi, and let [Image: pcbi.1005856.e276.jpg]D1S denote an [Image: pcbi.1005856.e277.jpg]IG(2bS/νS*,2bS/νS*3)-distributed random stable dominance time. Further, let [Image: pcbi.1005856.e278.jpg]pSS* denote the probability to stay in the stable state (8), and let Yi denote the hidden state during the i-th dominance time. Then it holds

Now we integrate across all possible values t that [Image: pcbi.1005856.e280.jpg]D1S can take, using independence of B and P during [Image: pcbi.1005856.e281.jpg]D1S, IG distribution of [Image: pcbi.1005856.e282.jpg]D1S and normal distribution of [Image: pcbi.1005856.e283.jpg]BD1S. Further we note that the probability in the denominator equals [Image: pcbi.1005856.e284.jpg]pSS*, which yields

[Image: pcbi.1005856.e285.jpg]μS+=(pSS*)-1∫0∞f2bS/νS*,2bS/νS*3IG(t)(1-ΦνB*t,t(b˜S))tdt.(15)

Analogously, we obtain

[Image: pcbi.1005856.e286.jpg]μS−=E[D1S|Y1=S,Y2=U]=(1−pSS*)−1∫0∞f2bS/νS*,2bS/νS*3IG(t)ΦνB*t,t(b˜S)tdt.(16)

Expected relative time spent in the stable state

In this section we derive the formula for the expected relative time [Image: pcbi.1005856.e287.jpg]φS* spent in the stable state in the HBMi analyzed in the section ‘Application of the HBM to the sample data set’. To that end, let [Image: pcbi.1005856.e288.jpg]Nj*, j ∈ {S, U}, denote the (random) number of dominance times of the HBMi in state j before a state change. As the number of dominance times before a state change is geometrically distributed with probability [Image: pcbi.1005856.e289.jpg]1-pjj* its expectation is given by [Image: pcbi.1005856.e290.jpg]E[Nj*]=(1-pjj*)-1. (Again, due to taking the mean of different drifts, this derivation is only a close approximation, but we omit approximation signs here for convenience.) Again, we derive [Image: pcbi.1005856.e291.jpg]φS* as

The length of a stable phase is a random variable distributed like [Image: pcbi.1005856.e293.jpg]∑i=1NS*DiS, where [Image: pcbi.1005856.e294.jpg]D1S,…,DNS*S denote the random stable dominance times. Given [Image: pcbi.1005856.e295.jpg]NS*, we use linearity of expectation to compute the conditional expectation

Now we take expectation over [Image: pcbi.1005856.e297.jpg]NS* to find

An analogous result holds for the expected length of an unstable phase. This yields

[Image: pcbi.1005856.e299.jpg]φS*=pSS*1-pSS*μS++μS-pSS*1-pSS*μS++μS-+pUU*1-pUU*μU++μU-.(17)

HBMi: Derivation of the alternation rate ρ as used in Fig 19B

Here, we derive a formula for the alternation rate ρ used in Fig 19 given the HBMi parameter set [Image: pcbi.1005856.e300.jpg](μS*,σS*,μU*,σU*,pSS*,pUU*). For each length of blank display lb and the value of νS, νU as shown in Fig 19A the mean drifts per second in the stable and the unstable state [Image: pcbi.1005856.e301.jpg]νS*, [Image: pcbi.1005856.e302.jpg]νU* and the mean drift of the background process [Image: pcbi.1005856.e303.jpg]νB* are derived using Eq (4). Then, we use Eqs (8) and (11) to derive the values [Image: pcbi.1005856.e304.jpg](μS*,σS*,μU*,σU*,pSS*,pUU*) given the mean drifts per second, and we recall Eq (17) for [Image: pcbi.1005856.e305.jpg]φS*, and analogously for [Image: pcbi.1005856.e306.jpg]φU*

We now show that the alternation rate can be described as

[Image: pcbi.1005856.e307.jpg]ρ≔limΔ→∞E[N*(Δ)]Δ=φS*μS*+φU*μU*,(18)

where N*(Δ) denotes the number of perceptual changes in an interval of length Δ. We split up [Image: pcbi.1005856.e308.jpg]N*(Δ)=NS*(Δ)+NU*(Δ), where [Image: pcbi.1005856.e309.jpg]NS*(Δ) and [Image: pcbi.1005856.e310.jpg]NU*(Δ) denote the number of perceptual changes in the respective stable and unstable phases in the interval of length Δ. We then show

and analogously for the unstable state. To that end, let ΔS be the time spent in the stable state in a time interval of length Δ. By the Elementary Renewal Theorem [e.g. 56] it holds

as Δ → ∞ naturally implies ΔS → ∞. According to the definition of [Image: pcbi.1005856.e313.jpg]φS* as the expected relative time spent in the stable state, we get [Image: pcbi.1005856.e314.jpg]ΔS/Δ→φS* in probability. This yields the claim

Experimental protocol and data preprocessing

Details on experimental protocol

The sample data set for main analysis was collected as described in [10], the sample data for the analysis of reproducibility was collected as described in [9]. In both cases, participants took part in a behavioral experiment including continuous and intermittent stimulation with the same ambiguous stimulus. This stimulus was a structure-from-motion stimulus that appears as a rotating sphere. The rotation direction of the sphere is ambiguous and equally compatible with leftward and rightward percepts. For continuous stimulation, the sphere was presented throughout four minutes in which participants’ perception spontaneously alternated between the two percepts. Participants indicated perceptual alternations via button presses on a keyboard. For intermittent stimulation, the sphere was presented for twenty minutes repeatedly for short intervals of 0.6 s interleaved by blank screens of 0.8 s duration. Here, participants indicated the perceived rotation direction at every 1.4 s at each presentation of the sphere. Please refer to [9, 10] for a detailed description of the data collection. As suggested in [10], missing responses during intermittent stimulation were replaced by their preceding responses because the reported percept typically persisted to the next available response.

Investigation of history dependence in the response patterns

Long-term dependencies in the data set of [10] during continuous presentation are analyzed using the Pearson correlation coefficient cH between the dominance times and the cumulative history H as introduced in [31]. The history H is a function of the length and recency of previously dominated percepts. For each of the 57 subjects with at least five dominance times, cH is estimated as explained in [31]. To assess statistical significance, 1000 data sets are obtained for each subject by permutation of the dominance times to approximate the distribution of cH under the null hypothesis of independent and identically distributed dominance times. Statistical significance on the 5% level is obtained by comparison of the empirical history cH to the 95% quantile of the distribution of cH derived from the permuted data sets.

Supplementary Material

Open In Web Browser

Acknowledgments

We thank Brooks Ferebee and Michael Messer for inspiring discussions.

Articles from PLoS Computational Biology are provided here courtesy of Public Library of Science

PMC Copyright Notice

The articles available from the PMC site are protected by copyright, even though access is free. Copyright is held by the respective authors or publishers who provide these articles to PMC. Users of PMC are responsible for complying with the terms and conditions defined by the copyright holder.

Users should assume that standard copyright protection applies to articles in PMC, unless an article contains an explicit license statement that gives a user additional reuse or redistribution rights. PMC does not allow automated/bulk downloading of articles that have standard copyright protection.

See the copyright notice on the PMC site, https://www.ncbi.nlm.nih.gov/pmc/about/copyright/, for further details and specific exceptions.

References

Necker LA. Observations on some remarkable optical phaenomena seen in Switzerland; and on an optical phaenomenon which occurs on viewing a figure of a crystal or geometical solid. London and Edinburgh Philosophical Magazine and Journal of Science. 1832;1:329–337.

von Helmholtz H. Treatise on Physiological Optics. vol. 3 Birmingham, AL: The Optical Society of America; 1866.

Levelt W. On binocular rivalry. Soesterberg, The Netherlands: Institute for Perception; 1965 Available from: https://pdfs.semanticscholar.org/7d66/4283485e4297924fe0abd03c7c75812c74a7.pdf.

Leopold DA, Logothetis NK. Multistable phenomena:changing views in perception. Trends Cogn Sci. 1999;3(7):254–264. doi: 10.1016/S1364-6613(99)01332-7 10377540

Hohwy J, Roepstorff A, Friston K. Predictive coding explains binocular rivalry: An epistemological review. Cognition. 2008;108:687–701. doi: 10.1016/j.cognition.2008.05.010 18649876

Sterzer P, Kleinschmidt A, Rees G. The neural bases of multistable perception. Trends Cogn Sci. 2009;13(7):310–318. doi: 10.1016/j.tics.2009.04.006 19540794

Braun J, Mattia M. Attractors and noise: Twin drivers of decisions and multistability. NeuroImage. 2010;52(3):740–751. doi: 10.1016/j.neuroimage.2009.12.126 20083212

Weilnhammer V, Stuke H, Hesselmann G, Sterzer P, Schmack K. A predictive coding account of bistable perception—a model-based fMRI study. PLOS Comput Biol. 2017;13(5):e1005536 https://doi.org/10.1371/journal.pcbi.1005536. 28505152

Schmack K, Gòmez-Carrillo de Castro A, Rothkirch M, Sekutowicz M, Rössler H, Haynes JD, et al Delusions and the Role of Beliefs in Perceptual Inference. J Neurosci. 2013;33(34):13701–13712. doi: 10.1523/JNEUROSCI.1778-13.2013 23966692

Schmack K, Schnack A, Priller J, Sterzer P. Perceptual instability in schizophrenia: Probing predictive coding accounts of delusions with ambiguous stimuli. Schizophr Res Cog. 2015;2(2):72–77. doi: 10.1016/j.scog.2015.03.005

Brascamp JW, Klink PC, Levelt WJ. The’laws’ of binocular rivalry: 50 years of Levelt’s propositions. Vision Res. 2015;109:20–37. doi: 10.1016/j.visres.2015.02.019 25749677

Levelt WJ. Note on the distribution of dominance times in binocular rivalry. Br J Psychol. 1967;58(1):143–145. doi: 10.1111/j.2044-8295.1967.tb01068.x 5582864

Moreno-Bote R, Rinzel J, Rubin N. Noise-induced alternations in an attractor network model of perceptual bistability. J Neurophysiol. 2007;98(3):1125–1139. doi: 10.1152/jn.00116.2007 17615138

Brascamp JW, Pearson J, Blake R, van den Berg AV. Intermittent ambiguous stimuli: Implicit memory causes periodic perceptual alternations. J Vis. 2009;9(3):3, 1–23. doi: 10.1167/9.3.3 19757942

Leopold DA, Wilke M, Maier A, Logothetis NK. Stable perception of visually ambiguous patterns. Nat Neurosci. 2002;5(6):605–609. doi: 10.1038/nn0602-851 11992115

Murata T, Matsui N, Miyauchi S, Kakita Y, Yanagida T. Discrete stochastic process underlying perceptual rivalry. NeuroReport. 2003;14(10):1347–1352. doi: 10.1097/00001756-200307180-00013 12876471

Wilson HR. Minimal physiological conditions for binocular rivalry and rivalry memory. Vision Res. 2007;47(21):2741–2750. doi: 10.1016/j.visres.2007.07.007 17764714

Gigante G, Mattia M, Braun J, Del Giudice P. Bistable Perception Modeled as Competing Stochastic Integrations at Two Levels. PLOS Comput Biol. 2009;5(7):e1000430 doi: 10.1371/journal.pcbi.1000430 19593372

Gershman SJ, Vul E, Tenenbaum JB. Multistability and Perceptual Inference. Neural Comput. 2012;24(1):1–24. doi: 10.1162/NECO_a_00226 22023198

Pastukhov A, Garcia-Rodriguez P, Haenicke J, Guillamon A, Deco G, Braun J. Multi-stable perception balances stability and sensitivity. Front Comput Neurosci. 2013;7(17). doi: 10.3389/fncom.2013.00017 23518509

Cao R, Pastukhov A, Mattia M, Braun J. Collective Activity of Many Bistable Assemblies Reproduces Characteristic Dynamics of Multistable Perception. J Neurosci. 2016;36(26):6957–6972. doi: 10.1523/JNEUROSCI.4626-15.2016 27358454

Orbach J, Ehrlich D, Heath H. Reversibility of the Necker cube. I. An examination of the concept of’satiation or orientation’. Percept Mot Skills. 1963;17:439–458. doi: 10.2466/pms.1963.17.2.439 14065532

Maier A, Wilke M, Logothetis NK, Leopold DA. Perception of temporally interleaved ambiguous patterns. Curr Biol. 2003;13(13):1076–1085. doi: 10.1016/S0960-9822(03)00414-7 12842006

Kornmeier J, Bach M. Early neural activity in Necker-cube reversal: Evidence for low-level processing of a gestalt phenomenon. Psychophysiology. 2004;41(1):1–8. doi: 10.1046/j.1469-8986.2003.00126.x 14692995

Pastukhov A, Braun J. Structure-from-motion: dissociating perception, neural persistence, and sensory memory of illusory depth and illusory rotation. Atten Percept Psychophys. 2013;75(2):322–340. doi: 10.3758/s13414-012-0390-0 23150214

Cao R, Braun J, Mattia M. Stochastic Accumulation by Cortical Columns May Explain the Scalar Property of Multistable Perception. Phys Rev Lett. 2014;113(9):098103 doi: 10.1103/PhysRevLett.113.098103 25216009

Laing CR, Chow CC. A Spiking Neuron Model for Binocular Rivalry. J Comput Neurosci. 2002;12(1):39–53. doi: 10.1023/A:1014942129705 11932559

Noest AJ, van Ee R, Nijs MM, van Wezel RJA. Percept-choice sequences driven by interrupted ambiguous stimuli: A low-level neural model. J Vis. 2007;7(8):10, 1–14. doi: 10.1167/7.8.10 17685817

Walker P. Stochastic properties of binocular rivalry alternations. Percept Psychophys. 1975;18(6):467–473. doi: 10.3758/BF03204122

Lehky SR. Binocular rivalry is not chaotic. Proc R Soc Lond B Biol Sci. 1995;259:71–76. doi: 10.1098/rspb.1995.0011

Pastukhov A, Braun J. Cumulative history quantifies the role of neural adaptation in multistable perception. J Vis. 2011;11(10):12, 1–10. doi: 10.1167/11.10.12 21931128

Seshadri V. The Inverse Gaussian Distribution. Oxford University Press; 1993.

Baum LE, Petrie T, Soules G, Weiss N. A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains. Annals of Mathematical Statistics. 1970;41(1):164–171. doi: 10.1214/aoms/1177697196

Rabiner LR. A tutorial on hidden markov models and selected applications in speech recognition. IEEE Trans Inf Theory. 1989;77(2):257–284.

Bilmes JA. A gentle tutorial of the EM algorithm and its application to parameter estimation for gaussian mixture and hidden markov models. International Computer Science Institute; 1998.

MacDonald IL, Zucchini W. Hidden Markov and other models for discrete-valued time series vol. 70 of Monographs on Statistics and Applied Probability. Chapman&Hall, London; 1997.

Ephraim Y, Merhav N. Hidden Markov processes. IEEE Trans Inf Theory. 2002;48(6):1518–1569. doi: 10.1109/TIT.2002.1003838

Samanta M. On Tests Of Equality Of Two Inverse Gaussian Distributions. S Afr Stat J. 1985;19(2):83–95.

Barndorff-Nielsen OE. Normal Inverse Gaussian Distributions and Stochastic Volatility Modeling. Scand Stat Theory Appl. 1997;24(1):1–13.

Turner R. Direct maximization of the likelihood of a hidden Markov model. Comput Stat Data Anal. 2008;52(9):4147–4160. doi: 10.1016/j.csda.2008.01.029

Stuke H, Stuke H, Weilnhammer VA, Schmack K. Psychotic Experiences and Overhasty Inferences Are Related to Maladaptive Learning. PLOS Comput Biol. 2017;13(1):e1005328 doi: 10.1371/journal.pcbi.1005328 28107344

Peters E, Joseph S, Garety P. Measurement of delusional ideation in the normal population: introducing the PDI (Peters et al. Delusions Inventory). Schizophr Bull. 1999;25(3):533–576. doi: 10.1093/oxfordjournals.schbul.a03340110478787

Fletcher PC, Frith CD. Perceiving is believing: a Bayesian approach to explaining the positive symptoms of schizophrenia. Nat Rev Neurosci. 2009;10:48–58. doi: 10.1038/nrn2536 19050712

Evans SL, Averbeck BB, Furl N. Jumping to conclusions in schizophrenia. Neuropsychiatr Dis Treat. 2015;11:1615–1624. doi: 10.2147/NDT.S56870 26170674

Ross RM, McKay R, Coltheart M, Langdon R. Jumping to Conclusions About the Beads Task? A Meta-analysis of Delusional Ideation and Data-Gathering. Schizophr Bull. 2015;41(5):1183–1191. https://doi.org/10.1093/schbul/sbu187. 25616503

Adams RA, Stephan KE, Brown HR, Frith CD, Friston KJ. The computational anatomy of psychosis. Front Psychiatry. 2013;4:47 doi: 10.3389/fpsyt.2013.00047 23750138

van Ee R. Dynamics of perceptual bi-stability for stereoscopic slant rivalry and a comparison with grating, house-face, and Necker cube rivalry. Vision Res. 2005;45(1):29–40. doi: 10.1016/j.visres.2004.07.039 15571736

Rankin J, Sussman E, Rinzel J. Neuromechanistic Model of Auditory Bistability. PLOS Comput Biol. 2015;11(11):e1004555 https://doi.org/10.1371/journal.pcbi.1004555. 26562507

Brascamp JW, van Ee R, Noest AJ, Jacobs RAH, van den Berg A. The time course of binocular rivalry reveals a fundamental role of noise. J Vis. 2006;6(11):1244–1256. https://doi.org/10.1167/6.11.8. 17209732

Pastukhov A, Braun J. Perceptual reversals need no prompting by attention. J Vis. 2007;7(10):5, 1–17. doi: 10.1167/7.10.5 17997674

Dieter KC, Brascamp J, D T, Blake R. Does visual attention drive the dynamics of bistable perception? Atten Percept Psychophys. 2016;78(7):1861–1873. doi: 10.3758/s13414-016-1143-2 27230785

Grimmett G, Stirzaker D. Probability and random processes. 3rd ed Oxford University Press, New York; 2001.

Wald A. On cumulative sums of random variables. The Annals of Mathematical Statistics. 1944;15(3):283–296. doi: 10.1214/aoms/1177731235

Dennis JE, Schnabel RB. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Prentice-Hall, Englewood Cliffs, NJ; 1983.

Powell MJD. A direct search optimization method that models the objective and constraint functions by linear interpolation In: Gomez S, Hennart JP, editors. Advances in Optimization and Numerical Analysis; 1994 p. 51–67.

Ross SM. Stochastic processes. Wiley, New York; 1996.

[Back]

Fig 1. Examples of response patterns to a bistable stimulus.

Response patterns to continuous (green, A,B) and intermittent (blue, C-F) presentation from the dataset reported in [10]. While the distribution of dominance times tends to be unimodal in the continuous case, stable and unstable phases interchange in intermittent stimulation. In addition, response patterns can be highly variable across subjects.

[Back]

[Back]

Fig 2. Alternation rates in control subjects and subjects with schizophrenia.

During continuous presentation healthy controls (C) showed higher alternation rates compared to patients with schizophrenia (S) (left), while the opposite could be observed for intermittent presentation (right). Each grey dot indicates the perceptual alternation rate from one individual participant, colored diamonds and lines indicate group medians with 25%/75% quantiles. Two-sided Wilcoxon tests yielded p < .1 for both continuous and intermittent stimulation.

[Back]

[Back]

Fig 3. A simple HMM for bistable perception.

(A) One state describes a unimodal distribution of dominance times under continuous presentation. (B) Two states (stable, S, and unstable, U) produce long and short dominance times under intermittent presentation.

[Back]

[Back]

Fig 4. Precision of parameter estimates in the one-state HMM.

For each of 61 parameter constellations in continuous presentation, 1000 simulations were performed with sample sizes as in the original data. (A) Median of the relative error (RE) for each parameter. (B) The CV as a function of the number n of simulated dominance times. Black points indicate constellations with mean RE across the parameters smaller than 0.25.

[Back]

[Back]

Fig 5. Precision of parameter estimates in the two-state HMM.

For each of 61 parameter constellations in intermittent presentation, 1000 simulations were performed. Log(median) of RE (for μS, σS, μU, σU) or of AE (for pSS, pUU) for T₁ = 1200 (A) and T₂ = 3600 (B). Black lines indicate constellations with mean errors across the parameters <0.25.

[Back]

[Back]

Fig 6. Comparison of empirical response patterns to patterns simulated with the HMM.

Examples of response patterns to continuous (green, panels A and B) and intermittent (blue, C-F) stimulation repeated from Fig 1 and corresponding simulations within the HMM (orange). The parameter estimates are given in Tables 1 and 2, respectively.

[Back]

[Back]

Table 1. Estimated HMM parameters of the response pattens to continuous presentation shown in Fig 1A and 1B.

com. | [Image: pcbi.1005856.e007]μ^ | [Image: pcbi.1005856.e008]σ^
A | 10.50 | 8.18
B | 6.69 | 3.58

DOI: 10.1371/journal.pcbi.1005856.t001

[Back]

[Back]

Table 2. Estimated HMM parameters of the response pattens to intermittent presentation shown in Fig 1C–1F.

com. | [Image: pcbi.1005856.e009]μ^S | [Image: pcbi.1005856.e010]σ^S | [Image: pcbi.1005856.e011]μ^U | [Image: pcbi.1005856.e012]σ^U | [Image: pcbi.1005856.e013]p^SS | [Image: pcbi.1005856.e014]p^UU
C | 186.45 | 30.50 | 5.01 | 3.06 | 0.67 | 0.96
D | 67.12 | 50.81 | 3.25 | 2.76 | 0.33 | 0.82
E | 372.28 | 68.30 | NA | NA | 1.00 | NA
F | 77.13 | 9.99 | 5.37 | 6.26 | 0.00 | 0.99

DOI: 10.1371/journal.pcbi.1005856.t002

[Back]

[Back]

Fig 7. Comparison of distribution of dominance times with theoretical distribution.

The theoretical IG distribution in the HMM fitted to the empirical distribution of dominance times for continuous (A) and intermittent (B) presentation shown in Fig 6A and 6C.

[Back]

[Back]

Fig 8. Reproducibility of response patterns.

Parameter estimates of log(μ) (A) and log(σ) (B) of the IG distribution in two sessions with the same individuals, dataset reported in [9]. The logarithm was applied due to asymmetric distributions of the parameter estimates. Stars indicate highly significant (p < .0001) correlation of parameter estimates across different sessions.

[Back]

[Back]

Fig 9. Differences in the HMM parameter estimates between subjects with schizophrenia and control subjects.

(A) [Image: pcbi.1005856.e028]μ^ during continuous presentation and [Image: pcbi.1005856.e029]p^SS and [Image: pcbi.1005856.e030]φ^S for intermittent presentation (B). Each grey dot indicates one individual participant’s data, colored diamonds indicate medians, horizontal bars indicate 25%/75%-quantiles.

[Back]

[Back]

Fig 10. The HBMc.

The first passage times (Hi) (green) of a Brownian motion Pt (black, panel A) with drift ±ν₀ at borders ±b indicate the times of the percept changes (orange, panel B). The Brownian motion is assumed to summarize the activity difference of two conflicting neuronal populations with only two parameters.

[Back]

[Back]

Fig 11. The perception process P in the HBMi during intermittent presentation.

During presentation, P has drift ν₀. During blank displays (yellow), P has drift νS in the stable phase (A), and drift νU in the unstable phase (B). Typically, we have νS ≤ ν₀ and νU ≥ ν₀. The borders are bS (light blue horizontal line) in the stable phase and bU (blue horizontal line) in the unstable phase.

[Back]

[Back]

Fig 12. The HBMi.

The perception process P, the background process B and the resulting percept. A: A simulation on [0, 500]. B: The same realization, zoomed in on the time interval [408, 458]. Stable phases indicated by white background, unstable phases indicated by gray background. The beginnings of stable and unstable phases are marked with light blue and blue arrows, respectively.

[Back]

[Back]

Fig 13. Motivation of HBMi model assumptions.

For explanation and details see text.

[Back]

[Back]

Fig 14. Impact of HBMi parameter values on the response patterns.

Examples of simulated response patterns are shown for different values of the three quantities [Image: pcbi.1005856.e071]CVS*, [Image: pcbi.1005856.e072]2bSνB*/b˜SνS* and [Image: pcbi.1005856.e073]b˜U. The quantities for panels A-H were [Image: pcbi.1005856.e074]CVS*={0.2,0.2,0.2,0.2,1,1,1,1}, [Image: pcbi.1005856.e075]2bSνB*/b˜SνS*={4,4,0.4,0.4,4,4,0.4,0.4}, [Image: pcbi.1005856.e076]b˜U={0,3,0,3,0,3,0,3}.

[Back]

[Back]

Fig 15. Precision of parameter estimates in the HBMi.

For each of the 61 parameter constellations estimated from the sample data set, 1000 simulations were performed with the HBMi. (A) and (B): log(median(REs)) of the border parameters bS, bU, [Image: pcbi.1005856.e144]b˜S, [Image: pcbi.1005856.e145]b˜U, and log(median(AEs)) of the drift parameters [Image: pcbi.1005856.e146]νS*, [Image: pcbi.1005856.e147]νU*, [Image: pcbi.1005856.e148]νB* with T = 1200s (A) and T = 3600s (B). (C) and (D) log(median(REs)) for the derived model parameters [Image: pcbi.1005856.e149]μS*, [Image: pcbi.1005856.e150]σS*, [Image: pcbi.1005856.e151]μU*, [Image: pcbi.1005856.e152]σU* and log(median(AEs)) for the parameters [Image: pcbi.1005856.e153]pSS* and [Image: pcbi.1005856.e154]pUU* for T = 1200s (C) and T = 3600s (D). Parameter combinations with mean errors < 0.25 across all variables plotted in black.

[Back]

[Back]

Fig 16. HBM: Examples of simulated response patterns.

Example empirical response patterns from Fig 1 and response patterns simulated by the HBM (orange). Responses to continuous and intermittent presentation are plotted in green and blue, respectively. The estimated parameters used for simulation are given in Tables 3 and 4, respectively. (A) and (B) Continuous presentation and HBMc. (C)-(F): Intermittent presentation and HBMi. Below the orange response patterns, one can see the perception process P and (for the HBMi) the background process B corresponding to the respective simulation.

[Back]

[Back]

Table 3. Estimated HBMi parameter combinations of the typical response pattens to continuous presentation shown in Fig 1A and 1B.

com. | [Image: pcbi.1005856.e155]b^ | [Image: pcbi.1005856.e156]ν^0
A | 2.08 | 0.40
B | 2.42 | 0.72

DOI: 10.1371/journal.pcbi.1005856.t003

[Back]

[Back]

Table 4. Estimated HBMi parameter combinations of the typical response pattens to intermittent presentation shown in Fig 1C–1F.

com. | [Image: pcbi.1005856.e157]ν^0 | [Image: pcbi.1005856.e158]b^S | [Image: pcbi.1005856.e159]ν^S | [Image: pcbi.1005856.e160]b^U | [Image: pcbi.1005856.e161]ν^U | [Image: pcbi.1005856.e162]b˜^S | [Image: pcbi.1005856.e163]b˜^U | [Image: pcbi.1005856.e164]ν^B
C | 0.72 | 41.74 | 0.24 | 1.84 | 0.74 | 49.71 | 2.39 | 0.54
D | 0.27 | 5.42 | 0.08 | 1.06 | 0.93 | 4.56 | 1.13 | 0.06
E | 0.32 | 52.58 | 0.25 | NA | NA | 0.00 | NA | 17.5
F | 0.57 | 13.5 | 0.24 | 1.00 | 0.24 | 110.28 | 0.77 | 1.90

DOI: 10.1371/journal.pcbi.1005856.t004

[Back]

[Back]

Fig 17. Mean dominance times in stable state as a function of the successive state.

In the empirical data set, the mean dominance time before a state transition SU is shorter than the mean dominance time before SS (A). This observation can be reproduced in the HBMi (C), while the expected dominance time in the HMM (B) is equally long for intervals with and without transition, i.e., independent from transition between states. In the empirical data set [10], the state transitions were estimated using the Viterbi paths [34]. Analogous results were obtained using a fixed threshold. For the HBMi, the expected dominance times were derived from the HBMi parameters according to Eq (15) in the Materials and methods. In the HMM, expected dominance times correspond to [Image: pcbi.1005856.e165]μ^S.

[Back]

[Back]

Fig 18. Group differences between patients with schizophrenia and control subjects in the HBM.

The HBMc parameter [Image: pcbi.1005856.e169]b^ (A) and the relation [Image: pcbi.1005856.e170]b^/b^S (B). (C): Differences in [Image: pcbi.1005856.e171]b^S/b˜^S and [Image: pcbi.1005856.e172]p^SS*. The raw data together with the median (colored diamonds) and 25%- and 75%-quantiles are shown. All p-values of a two-sided Wilcoxon test were below.1.

[Back]

[Back]

Fig 19. HBMi extension for different lengths of blank display.

Parameters derived for subject B and C in Tables 3 and 4, which is the same subject during continuous and intermittent presentation. A) The mean drifts νS, νU decrease with the length of the blank display lb. B) The resulting mean alternation rate ρ under intermittent stimulation (black) is derived using formula (18) (Materials and methods). The empirical mean alternation rates per minute of the subject during intermittent viewing with lb = 0.8 and during continuous viewing are marked by blue and green lines, respectively.

[Back]