A hierarchical stochastic model for bistable perception

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.


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].
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 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.
presentation, the potential neuronal mechanisms underlying the transitions between the different response properties remain unclear.
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.

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 d i , i = 1, 2, . . ., n. In the unimodal continuous case, we assume the d i to be the realizations of independent and identically (i.i.d.) distributed random variables D i , 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.
Formally, let Y ≔ (Y i ) i = 1,. . .,n describe a Markov chain on {S, U}, where S and U denote the stable and unstable state, respectively. Let p SU = 1 − p SS and p US = 1 − p UU denote the transition probabilities. The dominance times (d i ) i2{1,2,. . .,n} are assumed to be Inverse Gaussian distributed and conditionally independent given Y, with mean and standard deviation given by 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 1 , . . ., d n , assuming these are   [32] with mean μ and standard deviation σ, the maximum likelihood (ML) estimators are given bym 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 ðm i ;ŝ i Þ i¼1;...;61 with the parameter values underlying the simulation using the relative error (RE), defined e.g. for μ as REðmÞ ≔ jm À mj=m. 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).
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 1 = 1200 s and T 2 = 3600 s. For p SS and p UU the absolute errors (i.e., AEðp SS Þ ≔ jp SS À p SS j) are presented due to the small values of the two parameters. For the time horizon of the data, T 1 (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 2 (panel B), almost all parameter combinations showed errors smaller than 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 " t ¼ 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.
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).
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). 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 1 and n 2 , i.e. H 0 : μ 1 = μ 2 and σ 1 = σ 2 . The likelihood ratio derives as Under H 0 the quantity Q Ã n ≔ À 2ð1 À 1=6½1=n 1 þ 1=n 2 À 1=½12nÞ log Q n is approximately chi-square distributed with two degress of freedom. Thus, the test rejects the null hypothesis at level 5% if Q Ã n exceeds the 95%-th-quantile of the χ 2 (2)-distribution. In the  Tables 1 and 2 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(μ 1 ) and log(μ 2 ) (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 timem ( 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, E½length of a stable phase E½length of a stable phase þ length of an unstable phase 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 probabilitŷ p SS to stay in the stable state was higher in healthy controls. In addition, the mean dominance timem 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 form,m U ,p SS andφ S , two-sided Wilcoxon test). In the following section, we present a model that links these observed differences to potential underlying neuronal mechanisms.

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  [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.
https://doi.org/10.1371/journal.pcbi.1005856.g008 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 ±ν 0 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 ν 0 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 > 0 be a drift and T > 0 a time horizon, and let (W t ) t2[0,T] be a standard Brownian motion. The perception process P ≔ (P t ) t2[0,T] is then defined by The perception at time t ! 0 takes the value L (left) if S t = 1 and R (right) if S t = −1 and switches at the first-hitting times (H i ) i of the borders ±b defined by H 0 ≔ 0 and 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 (H i ) i , which also marks switches in the percept (panel B).
Parameter estimation in the HBMc makes use of the fact that the resulting dominance times are independent and IG distributed, with a known relation between the HBMc parameters b and ν 0 and the parameters μ and σ of the IG distribution, given by μ = 2b/ν 0 and s ¼ ffiffiffiffiffiffiffiffiffiffiffi 32]. Note that the CV of the dominance time distribution is therefore given by 1/(2bν 0 ). Thus, an increase in the border b also increases the mean dominance time and decreases the CV. An increase in the drift ν 0 decreases the mean dominance time, while again decreasing the CV.
The ML estimatorsb andn 0 can be derived viâ wherem andŝ 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 ν 0 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 b U at which perception and drift direction change is assumed smaller than the border b S 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.
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 ν 0 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 l p and blank display of length l b , PR and BL are given by The perception process (P t ) t is then given by whereỸ t 2 fS; Ug denotes the hidden state and (W t ) t denotes a standard Brownian motion. As a result, the mean drift per second is given by for states S and U, respectively. Because the periods l b and l p are typically short in relation to a dominance time, the behavior of P can be approximated by a Brownian motion with absolute drifts n Ã S and n Ã U , respectively. As in the HBMc, the sign of the drift S t ≔ S(P t , t) changes at We initialize P 0 ¼ À bỸ 0 for the initial stateỸ 0 which is the stable state with probability p S ≔ PðỸ 0 ¼ SÞ. The perception then takes the value L if S t = 1 and R if S t = −1 and switches at the first-hitting times (H i ) i of the borders ±b i comparable to Eq (2). Note that perception also changes during blank display. The dominance times are therefore again given by 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 ≔ (B t ) 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., where ðW t Þ t is a Brownian motion independent of (W t ) t . Again, the mean drift across PR and BL intervals is n Ã B ≔ l b Án B l b þl p . 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, b U < b S determine this switching as follows (see Fig 12).
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 ±b S in the stable state and between ±b U in the unstable state, the value of P is reset when the state changes, to the value sgn(P)b S when changing to the stable state and to sgn(P)b U 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. 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, b S > b U for the stable and the unstable state. This suggests different population sizes of neurons involved in the stable and unstable state. Typically b S > b > b U , 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,b S andb 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 (aboveb S ), the activity of S can increase rapidly again, causing another stable dominance time. In contrast, in case of weak previous activation (belowb 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 borderb 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 ðb S ; n Ã S Þ and ðb U ; n Ã U Þ can be interpreted analogously to the parameters (b, ν 0 ) in continuous stimulation. That means, an increase in the border (b S or b U ) increases the mean dominance time and decreases the CV in the respective state. An increase in the drift (n Ã S or n Ã 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 ; respectively: illustrates examples with small CV Ã S (panels A-D) and large CV Ã S (panels E-H). Second, the parametersb S and n Ã B can be interpreted best when compared to b S and n Ã S as follows. Consider the expected fraction ofb S reached by the background process at the end of a which is related to the transition probability from stable to unstable state. In case of a small background borderb S < b S and small n Ã S , the probability of B crossingb 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 parametersb U and n Ã U . Third, the parameterb 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 ofb U implies a low probability to reachb 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 ofb 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 ðb S ; n Ã 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 n Ã S (or n Ã U ), instead of exactly assuming drift ν 0 during stimulation and ν S (or ν U ) during blank display. However, the marginal distribution of P at multiples of such intervals l p + l b is identical to the marginal distribution of a Brownian motion with drift n Ã S (or n Ã U ) at these time points, and the differences can only be observed in the meantime. Because dominance times usually span multiple trials of duration l p + l b , the approximation is very close. As another similarity, both models use additional parameters ((p SS , p UU ) and ðb S ;b U ; n 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 − p SS ) and (1 − p UU ) 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 reachedb S at the end of the respective dominance time. Therefore, the transition probabilityp SU ðd i Þ depends on the duration d i 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 d i is given by an increment of a Brownian motion with drift n Ã B in the fixed time interval d i . Therefore the position is normally distributed with mean d i Á n Ã B and variance d i and the probability to remain in the stable state (which is the probability that the background process exceedsb S ) is given bỹ where F μ,σ 2() denotes the distribution function of the normal distribution with mean μ and variance σ 2 and Y i is the hidden state of the i-th dominance time. Analogously, for the transition from unstable to stable state, the transition probabilityp US ðd i Þ tends to decrease with longer dominance times, and we find Note that we use the approximate sign '%' because the drift of B is not exactly n Ã B throughout, but is assumed to change between ν 0 and ν B during stimulation and blank display, respectively, yielding a mean drift of n Ã 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 p SS and p UU in the HMM, we derive the marginal transition probabilities in the HBMi as the expected value ofp SS andp UU . As shown by [39], the positions X S and X U 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 where X S is NIG-distributed with parameters (0, . One should note that due to the difference in transition probabilities of the two models, the parameters ðb S ; n Ã 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 where the forward variable α j (i) denotes the probability of being in state j at time i while observing (d 1 , . . ., d i ). 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 with π j as starting distribution and f S and f U denoting the densities of IGð2b S =n Ã S ; Þ distributed random variables, respectively. Details on the maximization algorithm can be found in the Materials and methods. After estimation of Θ, the estimates ofn S andn U can be obtained using Eq (4), and analogously forn 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 b S , b U ,b S ,b U , we use the RE, while for the typically small drift parameters n Ã S , n Ã U , n Ã B , the AE is used. Again, one set of simulations uses the time horizon of the empirical data, T 1 = 1200s (Fig 15A), and a second simulation was performed using T 2 = 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 ðb S ; b U ;b S ;b U ; n Ã S ; n Ã U ; n Ã 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 1 and still in only 44 out of 61 for the tripled sample size of T 2 . 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 ðm Ã S ; s Ã S Þ and ðm Ã U ; s Ã 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 and analogously for m Ã U and s Ã U , where the approximation is again due to the minimal difference between the mean drift n Ã S and the changing drift ν S + ν 0 during presentation and blank display. In addition, we consider the transition probabilities between the states, p Ã SS and p Ã UU (cmp. Eq (8)). Fig 15C and 15D show the median REs for m Ã S , s Ã S , m Ã U , s Ã U and median AEs for p Ã SS , p Ã UU obtained in the 1000 simulations of length T 1 (panel C) and T 2 (panel D). Concerning this parameterization, 51 parameter combinations yielded average errors smaller than 0.25 across all variables for T 1 , while as many as 58 out of 61 cases showed errors smaller than 0.25 for T 2 . These simulations suggest that the parameterization ðm Ã S ; s Ã S ; m Ã U ; s Ã U ; p Ã SS ; p Ã UU Þ yields more reliable parameter estimates in the HBMi than the original parameterization of borders and drifts.
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).
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).
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 timem in the HMM (Fig 9A) therefore equals the corresponding value 2b=n 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 2b=n 0 . When analyzing the individual parameters b and ν 0 , we found no group differences in the drift ν 0 , 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).
Second, regarding intermittent presentation, we focussed on the derived parameters ðm Ã S ; s Ã S ; m Ã U ; s Ã U ; p Ã SS ; p Ã UU Þ, 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, φ Ã S (p < .1, two-sided Wilcoxon test, see Eq (17) in the Materials and methods). Again, the probability p Ã SS 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 p Ã SS . Noting that no difference was observed in the mean stable dominance time m Ã S ¼ 2b S =n Ã S and in the relation n Ã S =n Ã B of the drift in stable state and the drift of the background process, one interesting parameter is b S =b S . Keeping all other parameters constant, an increase in this quantity means an increase in p Ã SS because the crossing ofb 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 b S =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 b S /b, which describes a transition of population sizes from continuous to intermittent stimulus processing. In the data set, we observed increased values of b S /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 b S =b S and b S /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, 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 tom S . https://doi.org/10.1371/journal.pcbi.1005856.g017 Table 3. Estimated HBMi parameter combinations of the typical response pattens to continuous presentation shown in Fig 1A and 1B is a particular advantage of the HBM over the HMM, because the HBM provides mechanistic explanations and variables with potential neurophysiological interpretations.

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  Table 4. Estimated HBMi parameter combinations of the typical response pattens to intermittent presentation shown in Fig 1C- 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 ν 0 [21] between two borders ±b, where the first hitting times of the borders indicate percept changes. Roughly speaking, the drift ν 0 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 b S in the stable state to fast drift ν U and small border b U 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 b S ) 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, b S , 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 (b S >> b U ), 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 longterm 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 l b l p Á ν 0 /ν S . For longer blank displays, the mean drift of P during stable states, n Ã 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 b S (or b U ) 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 b S 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 l p 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 l b , 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 l b . 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 ν 0 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 l p (or 'ON'-period) and the length of the blank display l b (or 'OFF'-period). Particularly, the dependence between l b 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 (l b ) and ν U (l b ) are functions of the length of the blank display l b . In Fig 19A these mean drifts ν S , ν U decrease with l b , 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 l b > 0.7. This increased stability is caused first by a small drift ν S < ν 0 in that range. Second, it is also caused by the fact that the time interval l b in which the background process B has positive drift is longer, leading to an increased probability to reachb S and thus, to stay in the stable state. Estimation of the functions ν S (l b ) and ν U (l b ) from a suitable data set with variable lengths of blank displays would be an interesting task.
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.

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 1 , d 2 , . . ., d n ) the set of dominance times modeled as realizations of random variables D = (D 1 , . . ., D n ).
Continuous presentation: ML estimation. The HMM for continuous presentation assumes that all dominance times d i are independent and Inverse Gaussian distributed with mean μ and standard deviation σ. The log-likelihood function is then given by log Lðdjm; sÞ ¼ 1 2 Setting the partial derivatives with respect to μ and σ to zero yields the estimates of Eq (12) [32]:m Intermittent presentation: Baum-Welch-Algorithm. The HMM for intermittent presentation uses the parameter set ϑ = (μ S , σ S , μ U , σ U , p SS , p UU ) 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 1 , d 2 , . . ., d i 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 d i+1 , d i+2 , . . ., d n and being in state j at time i. To avoid underflow we normalize both variables [e.g. 34], resulting in the normalized variablesã j ðiÞ,b j ðiÞ. The normalized variables can be derived iteratively as follows where π j denotes the probability that the Markov chain starts in state j, p SU = 1 − p SS , p US = 1 − p UU and f IG m;s ðxÞ denotes the density of the IG distribution with expectation μ and standard deviation σ evaluated at x. Note that we suppress the dependence ofã j ðiÞ andb 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 1 , . . ., d n ) and the parameters ϑ g j ðijWÞ ¼ã j ðiÞb j ðiÞ a S ðiÞb S ðiÞ þã U ðiÞb U ðiÞ : 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 Y denote the set of all possible state sequences of the hidden Markov chain. Let Y = (Y 1 , . . ., Y n ) denote a Y-valued random variable and y = (y 1 , . . ., y n ) a realization of Y. i.e., the expectation of the complete-data log-likelihood L across all possible paths y 2 Y. The updated parameter set ϑ (m+1) is chosen such that it maximizes Q. For a fixed state sequence y = (y 1 , . . ., y n ) 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 PðY i ¼ y i jd; W ðmÞ Þ ¼ g y i ðijW m Þ and 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 In order to reduce the probability that the Baum-Welch-Algorithm will be captured in a local extremum, we chose ten equidistant values for m ðsÞ S ranging between 60 and 0.95max i d i , and for each value of m ðsÞ S we choose ten equidistant values for s ðsÞ S between 10 and 1:1m ð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 p SS ≔ 1. If the dominance time are only smaller than 30 seconds, we only estimate μ U and σ U by ML and use p UU = 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ŝ 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 withm S < 0:98m 15 wherem 15 ≔ ð1=jd1 d>15 jÞ P n i¼1 d i 1 d i >15 . This prevents dominance times smaller than 15 seconds to be considered for the estimation of μ S . Third, we requirem S < 1:02m 75 wherê m 75 ≔ ð1=jd1 d>75 jÞ P n i¼1 d i 1 d i >75 if any dominance time is larger than 75 seconds and m 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 N j , j 2 {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 − p jj , we have for the expectation Moreover, let ðD S i Þ i!1 be a sequence of independent IG(μ S , σ S )-distributed random variables and ðD U i Þ i!1 be a sequence of IG(μ U , σ U )-distributed random variables. We derive φ S as φ S ≔ E½length of a stable phase E½length of a stable phase þ length of an unstable phase : The length of a stable phase is a random variable distributed like . . . ; D S N S are independent from each other and also independent of N S . Therefore we have [53] E and analogously for the expected length of an unstable phase. This yields 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 ±ν 0 , where the sign of the drift changes at the first hitting time of the borders ±b. The two parameters (b, ν 0 ) 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/ν 0 and variance s 2 ¼ 2b=n 3 0 , or conversely, b ¼ ð1=2Þ Á p . This yieldŝ HBMi. The parameters of the HBMi are Y ¼ ðb S ; n Ã S ; b U ; n Ã U ;b S ;b U ; n Ã B Þ where ðb S ; n Ã S Þ and ðb U ; n Ã U Þ denote the border and drift parameters of the perception process P in the stable and the unstable state, respectively.b S ,b U , n Ã 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 LðdjYÞ % a S ðnÞ þ a 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 2 {S, U}, The likelihood then derives as Parameter estimation is then obtained by maximizing ∑ i log(c i ), 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 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 AÞ n U ! n S ; 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 b S , n Ã S are estimated like in the section 'HBMc' and b S is set to zero. In cases in which only dominance times up to 30 seconds are observed, we proceed analogously, where b U , n Ã U are estimated like in the section 'HBMc' and setb U ¼ 10 10 . In either case we set n Ã 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 m À S of a stable dominance time before a state change and its corresponding expected length m þ S before another stable dominance time. Let Y ¼ ðb S ; n Ã S ; b U ; n Ã U ;b S ;b U ; n Ã B Þ be the parameter set of an HBMi, and let D S 1 denote an IGð2b S =n Ã S ; ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2b S =n Ã S 3 p Þ-distributed random stable dominance time. Further, let p Ã SS denote the probability to stay in the stable state (8), and let Y i denote the hidden state during the i-th dominance time. Then it holds Now we integrate across all possible values t that D S 1 can take, using independence of B and P during D S 1 , IG distribution of D S 1 and normal distribution of B D S 1 . Further we note that the probability in the denominator equals p Ã SS , which yields Analogously, we obtain Expected relative time spent in the stable state. In this section we derive the formula for the expected relative time φ Ã 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 N Ã j , j 2 {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 1 À p Ã jj its expectation is given by E½N Ã j ¼ ð1 À p Ã jj Þ À 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 φ Ã S as E½length of a stable phase E½length of a stable phase þ length of an unstable phase : The length of a stable phase is a random variable distributed like P N Ã S i¼1 D S i , where D S 1 ; . . . ; D S N Ã S denote the random stable dominance times. Given N Ã S , we use linearity of expectation to compute the conditional expectation Now we take expectation over N Ã S to find An analogous result holds for the expected length of an unstable phase. This yields 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 ðm Ã S ; s Ã S ; m Ã U ; s Ã U ; p Ã SS ; p Ã UU Þ. For each length of blank display l b and the value of ν S , ν U as shown in Fig 19A the mean drifts per second in the stable and the unstable state n Ã S , n Ã U and the mean drift of the background process n Ã B are derived using Eq (4). Then, we use Eqs (8) and (11) to derive the values ðm Ã S ; s Ã S ; m Ã U ; s Ã U ; p Ã SS ; p Ã UU Þ given the mean drifts per second, and we recall Eq (17) for φ Ã S , and analogously for φ Ã

U
We now show that the alternation rate can be described as where N Ã (Δ) denotes the number of perceptual changes in an interval of length Δ. We split up 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 c H 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, c H 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 c H 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 c H to the 95% quantile of the distribution of c H derived from the permuted data sets.
Supporting information S1 Data. Raw data including dominance times from the studies [9,10]. Data include dominance times for intermittent (AlbertetalDataINT2015.RData) and continuous (AlbertetalData-CON2015.RData) stimulation for each of the 61 subjects of [10], and dominance times of two trials of continuous stimulation for each of the 105 subjects of [9] (AlbertetalDataCON2013. RData). A pdf file (AlbertetalDataDescription.pdf) provides a detailed description of the datasets. The data files and pdf file are packed in the supporting information file S1_Data.zip. (ZIP)