Table of Contents Author Guidelines Submit a Manuscript
Volume 2017 (2017), Article ID 4391587, 14 pages
Research Article

Impact of Time Delay in Perceptual Decision-Making: Neuronal Population Modeling Approach

1University of Warsaw, Stefana Banacha 2, 02-097 Warsaw, Poland
2Donders Institute for Brain, Cognition and Behavior, Kapittelweg 29, 6525 EN Nijmegen, Netherlands
3Radboud University Nijmegen Medical Centre, Geert Grooteplein Zuid 10, 6525 GA Nijmegen, Netherlands
4Nalecz Institute of Biocybernetics and Biomedical Engineering, Polish Academy of Sciences, Ks. Trojdena 4, 02-109 Warsaw, Poland

Correspondence should be addressed to Natalia Z. Bielczyk

Received 26 May 2017; Revised 15 July 2017; Accepted 25 July 2017; Published 6 September 2017

Academic Editor: Fathalla A. Rihan

Copyright © 2017 Urszula Foryś et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Impairments in decision-making are frequently observed in neurodegenerative diseases, but the mechanisms underlying such pathologies remain elusive. In this work, we study, on the basis of novel time-delayed neuronal population model, if the delay in self-inhibition terms can explain those impairments. Analysis of proposed system reveals that there can be up to three positive steady states, with the one having the lowest neuronal activity being always locally stable in nondelayed case. We show, however, that this steady state becomes unstable above a critical delay value for which, in certain parameter ranges, a subcritical Hopf bifurcation occurs. We then apply psychometric function to translate model-predicted ring rates into probabilities that a decision is being made. Using numerical simulations, we demonstrate that for small synaptic delays the decision-making process depends directly on the strength of supplied stimulus and the system correctly identifies to which population the stimulus was applied. However, for delays above the Hopf bifurcation threshold we observe complex impairments in the decision-making process; that is, increasing the strength of the stimulus may lead to the change in the neuronal decision into a wrong one. Furthermore, above critical delay threshold, the system exhibits ambiguity in the decision-making.

1. Introduction

Gamma-Aminobutyric Acid (GABA) is the most prevalent inhibitory neurotransmitter in the human brain [1, 2]. There is a body of evidence that aging has a major influence on the effectiveness of the GABAergic synapses [3]. Moreover, as found in the recent studies with the use of magnetic resonance spectroscopy (MRS, [4, 5]), local GABA concentrations in the frontal areas influence cognitive performance in aging adults. Taken together, aging can cause a reduced release of GABA to the intersynaptic cleft and decrease the quality of the synaptic transmission. This can result in an increase of the synaptic delays in the local inhibition.

On the other hand, in the aging process, the sensory capacity is declining, which affects the cognitive functions [68]. In particular, the working memory—which involves active manipulation of information—is affected, and this effect can further influence the perceptual decision-making [9]. (One important note here is that there is a difference between perceptual decision-making and making the abstract complex choices: the latter was reported not to be impaired in the elderly subjects [10], and some studies even report that elderly subjects are actually more efficient at making complex decisions [11, 12]. Therefore, in the paper we refer to the perceptual decision-making.).

In 1996, Salthouse [13] proposed a processing-speed theory of age-related deficits in cognition, for example, in working memory and perceptual decision-making. According to this theory, reduction in the processing speed can cause impairments in the cognitive functions for two major reasons. Cognitive functions can decline either because necessary logic operations cannot be executed within the time limit, or because the higher-order operations are blocked by slow execution of lower-order operations. This subject-matter was further experimentally investigated from different angles [1416], but the consensus in the field is that, in general, the elderly subjects are slower in perceptual decision-making than youngsters.

In this work, we propose a mechanism linking the two aspects of aging in cortical networks: the neurodegeneration in the local inhibitory synapses and the processing-speed related impairments in perceptual decision-making. This mechanism is based on a neuronal population model of decision-making based on a winner-take-all mechanism. The novelty lies in combining a winner-take-all mechanism well routed in the decision-making neuroscience, with the system of delayed differential equations representing the local inhibition within the two competing populations. With the use of this model, we are able to demonstrate that, for small synaptic delays in the local inhibition within the competing populations, the decision-making process depends directly on the strength of the stimulus, and the network is able to correctly identify the direction the stimulus came from. However, large delays can lead to a subcritical Hopf bifurcation resulting in complex decision-making process impairments. In particular, we demonstrate that, above the Hopf bifurcation point, increasing the strength of the stimulus can confuse the network and cause a wrong decision to be made. Furthermore, for delay values above critical threshold the system exhibits ambiguity in the decision-making. This effect can explain how the experimentally found difficulties in decision-making in elderly adults can be caused by loss in cognitive capacities [14].

The paper is organized in the following way. In Section 2, we introduce the model. In Section 3, qualitative analysis of the model, focusing on the stability analysis and Hopf bifurcation appearance, is made. Firstly, we present the analytic results for a symmetric model with no delay (Section 3.1) and then with positive delay (Section 3.2). In Section 4, we present the simulations undertaken to explore the dynamic repertoire of the model. In Section 5, we critically discuss the results and give recommendations for the future research.

2. Perceptual Decision-Making Model

The famous perceptual experiment on rhesus monkeys [17] by Shadlen and Newsome, involved a binary classification task: the monkeys had to assess whether the majority of the moving dots on the screen moved to the left or to the right. In the literature, one model proposed to model the decision-making in this experiment was the slow reverberation mechanism by Wang [18]. In Wang’s model, two populations of densely interconnected, spiking neurons compete with each other once being supplied by noisy inputs. This model involves simulations of two competing pools of stochastic spiking neurons. Since Wang’s work was published, some rodent [19] and computational [20] models were proposed to study perceptual choices.

In this work, we focus on modeling the most basic perceptual decision-making in neuronal networks. In such conditions, the network needs to disambiguate between two sensory stimuli. We consider changes in the firing rates of two positively interconnected self-inhibiting neuronal populations that receive external inputs and , respectively (Figure 1).

Figure 1: A neuronal population model of decision-making. Neurons in the first population project to the neurons in the second population and vice versa. Both neuronal populations receive self-inhibition with delay (, dashed arrows with flat heads). In addition, they receive external inputs and , respectively. The dynamics of this system is described by (1).

In order to describe the temporal dynamics of considered system (Figure 1), we introduce a system of delay differential equations (DDEs [21]) in the following form:where , denotes a model time scale, denotes the delay in self-inhibition, and is a coefficient describing the maximal capacity of a synapse, related to its anatomy. The function characterizes interactions between populations through synaptic plasticity which undergoes a Hebbian rule [22]; that is, the dynamics of synaptic weights is firing rate-dependent. For in-depth analytical and numerical investigations we have chosen a biologically plausible [23] sigmoid function .

In reality both populations considered in the proposed model are embedded in a larger network; therefore even in the absence of the population-specific sensory stimulus, they receive a constant input. Therefore, in the resting state, this system receives equal constant inputs to both nodes and, due to the system’s symmetry, both populations will be firing with the same rates. The symmetry breaks down if one of the populations receives an additional, external stimulation. The decision-making in this system means decoding which of the two populations received an additional stimulus (without estimation of the stimulus magnitude). The decoding is based on the difference between the two firing rates: the larger the difference , the more likely the decision that the population received the stimulation. The evidence behind each of the two options accumulates over time; therefore the psychometric function for the first population takes the integral form ofwhere is a parameter influencing the slope of the sigmoid function with respect to the cumulative difference . Similarly, we define the psychometric function for the second population asValues and can only asymptotically approach ; therefore we add a condition that if, at a given time point , exceeds the threshold value of (where is a given precision), the decision is made.

3. Qualitative Behavior of the Model

In this section, we provide the analytical results regarding the behavior of solutions of (1) for and constant symmetric input . Notice that the qualitative dynamics of (1) does not depend on .

3.1. Behavior of the Model for and Constant Symmetric Input

In this subsection, we present a detailed analysis of the model dynamics for and , as it is a crucial first step in the analysis of time-delayed models. While looking for steady states of (1), we need to solve the system of equationswhich yields . It is then obvious that both coordinates of any steady state are the same. Let us denote a steady state by . Clearly, , and the number of steady states depends on the shape of the graph of . Notice that the reference case is specific, as for the function is asymptotically linear , while for it tends to as .

Let us consider . Then we have , and therefore is increasing for , achieves its maximum at , and decreases to for . This implies that there are two steady states for and no steady state for , while for there is a bifurcation. The reference value , so there exist two steady states, and .

For , we have , and looking for zeros of we obtain , and we see that this quadratic equation has two positive solutions for , no real solution for , and one positive solution for . This means that for the function is increasing, for it has one maximum and tends either to (for ) or to (for , while for it is first increasing, then decreasing, and eventually increasing linearly to .

Corollary 1. For there is one steady state of (1); for there are two steady states for small values of and no steady state for larger values; for there is one steady state for small and sufficiently large values of , while for intermediate values 3 steady states exist.

Notice that due to its symmetric structure the system described by (1) always has solutions lying within a straight line . Clearly, assuming from both equations of (1) we obtainThe number of steady states determines the dynamics of (5). Let us assume that . Then there is only one steady state and the right-hand side of (5) is positive for and negative for . Therefore, is globally attractive. If , then there are two steady states for small , and the right-hand side is positive for and and negative for . This means that is locally stable, while is unstable. Moreover, solutions for tends to . If there is no steady state, then all solutions tend to , as they are increasing and unbounded. For there can be up to three steady states. Assume that there are three steady states . In this case and are stable, while is unstable. Moreover, for solutions tend to .

It is obvious that the symmetry of (1) implies that the phase space is divided into two symmetric subspaces by the straight line . Moreover, the dynamics of (5) is crucial in determining the whole model dynamics. Notice that and for we have . Similarly, for we have . Moreover, these straight lines are asymptotes for null-clines of (1). Therefore, asymptotic dynamics for large values of , is related to the dynamics of (5) for large , as it determines the direction of the vector-field in the region between the null-clines. Three generic types of the model dynamics when at least one steady state exists are presented in Figure 2.

Figure 2: Examples of the phase space portrait of (1) for and (a) (there is only one steady state, stable node at around ); (b) (there are three steady states: stable nodes at around and and a saddle at around ); (c) (there are two steady states: a stable node at around and a saddle at around ).
3.2. Model Behavior for and Constant Symmetric Input

Now, we aim to study the influence of the delay on the dynamics of the model. Our main goal is to show that there exists such time delay for which the steady state loses stability and a Hopf bifurcation occurs. Moreover, we would like to analyze the type of this bifurcation.

Before starting the analysis, we first state some general results that will be useful in this section. Let us consider a general DDEwith a smooth function . Let denote a characteristic function for (6) at . Assume thatwhere and are polynomials, . Together with (6) we considerfor which is a characteristic function.

Lemma 2. Assume that and have no common imaginary root and has a pair of purely imaginary simple eigenvalues for some critical value , and moreover these eigenvalues satisfy transversality condition . If , then is a pair of simple eigenvalues of for which satisfies transversality condition. Moreover, the eigenvalues of (6) cross imaginary axis in the same direction as the eigenvalues of (8).

Proof. We only need to check transversality condition. The derivative for (6) is obtained from the relationAs while , this implies which is the relation determining transversality condition for (8).

Lemma 2 holds under weaker assumptions, that is, for and being analytic functions, not only polynomials; compare [24] for more details. However, for our analysis, it is sufficient to consider the presented version.

The next lemma is a simple consequence of Proposition  1 from [24].

Lemma 3. Let , . Then there exists a pair of purely imaginary eigenvalues , , for   for which eigenvalues cross imaginary axis from left to right. Moreover, a steady state loses stability for and cannot gain it again for .

Proof. Following [24] we define an auxiliary function . It is obvious that is a simple zero of . Looking for critical delay related to this pair of eigenvalues , we need to solve a system of equationsIt is obvious that we obtain a sequence of critical delays , . However, as the derivative of the auxiliary function is positive, eigenvalues always cross imaginary axis from left to right (independently of ), which means that a switch of stability appears only at .

Now, we turn to the main topic of this subsection that is analysis of (1) for . While studying a Hopf bifurcation, we follow the ideas introduced in [25]. Let us rewrite the system described by (1) in its functional form (cf. [26, 27]); that is,where and for , is a linear part, and is nonlinear part of the system. Hence,where , , and . With the linear operator we are able to associate an operator which is a solution to (1) for initial data . In this way we obtain a strongly continuous semigroup generated by an infinitesimal generator (cf. [26, 27]). If for some critical value the generator has an eigenvalue , then a Hopf bifurcation occurs under the assumptions that the eigenvalues are simple and cross imaginary axis with nonzero speed when crosses .

Let us denote a characteristic matrix by ; that is,Looking for eigenvalues, we need to find zeros of the characteristic functionIt is obvious thatClearly, we can use Lemmas 2 and 3 in the analysis of stability switches for the steady state . In the considered case, the steady state is stable for ; that is, has negative zeros, yielding . Hence, both quasi-polynomials and satisfy assumptions of Lemma 3. This means that there are two sequences of critical delays and associated with and , respectively. However, the switch of stability can occur only for or , depending on the magnitude of those delays. Clearly, according to Lemma 2 eigenvalues cross imaginary axis in the same direction for both and , which means that they cross from left to right, and therefore the steady state loses stability for the smallest critical delay.

As a result, we can state the following theorem.

Theorem 4. Let . The steady state of (1) is locally asymptotically stable for and unstable for , and at a Hopf bifurcation occurs.

Proof. We only need to check that . Notice that . Moreover, the function , is decreasing. As we obtain

For the reference values of parameters for we obtain and . Considering we obtain and . For we have and

Hence, we study the bifurcation at and corresponding . Let us denote to shorten the notation. Therefore, . Notice that from the form of we obtain .

In the following, we base on the ideas presented in [25]. We know that is a purely imaginary eigenvalue of the infinitesimal generator for if there exists a vector such that and then is an eigenvector for at . Moreover, is also an eigenvalue for the adjoint operator for and is an eigenvector, where . In the considered case, we are able to choose such that , where is the derivative with respect to the first variable, here .

When looking for we obtainand thereforeClearly, and we can choose . Hence, is the eigenvector for the eigenvalue .

Next, we need to find a vector that satisfieswhich means that coordinates of satisfy the same relation as for , that is, . Next, calculatingwe obtain , that is, , and eventually

The type of the studied bifurcation is determined by a coefficient of the third term in Taylor expansion of an orbit. This coefficient readswhere denotes a first derivative with respect to the second variable , while , andwhere , denotes th derivative with respect to the first variable, (in fact is a constant function as it does not depend on ), and

If , then the bifurcation is supercritical; that is, periodic solutions appear for and are stable in such a case. If , then the bifurcation is subcritical: periodic solutions exist for , and as the steady state is stable in this case, periodic orbits are necessarily unstable.

First, we calculate the denominator of . We haveand we easily check that the real part of this expression is equal to .

Next, we would like to calculate the numerator of . To this end, we need to calculate the derivatives of the nonlinear part necessary to calculate , and we omit bar in the notation . Denoting by , , and test functions from , we obtainEvaluating the second and third derivative at we obtainwhere , is denoted by , to shorten the notation.

Using the formula above, we calculateTo calculate we need to evaluate . We haveMoreover,and thereforeFinally,To obtain the last parameter we first need to calculate .

Let us denote and calculate the first term of this matrix; that is,Hence,Next we calculate , whereFinally, we evaluateand therefore

Now, we are in a position to check the sign of , where . Clearly,and thereforewhere and .

Notice that, for small , the first term of the expression above is dominating, and therefore in the range of parameters we are interested in the fact that, as , the first term is positive and of the order of units, suggesting that . This implies that the bifurcation we study is subcritical, which therefore yields instability of appearing periodic orbits. However, in general it is necessary to calculate the exact value of , as it is difficult to guess its sign.

Now, we calculate for the reference parameters with . Let us recall that , , , , , and moreover . We obtain and . Consecutive terms in the brackets are equal to , , and , while the last fraction is . Eventually,which means that, in our reference case, the bifurcation is subcritical.

In order to complement the analysis above, we can also show that in the reference case the second bifurcation appearing for is subcritical as well. We again denote by , to shorten the notation. Now, the pair , satisfies yielding , and from this relation we obtain the vector . Hence, is the eigenvector for the eigenvalue . Next, we find such that . As before, coordinates of satisfies the same relation as for , that is, , . Calculating the denominator of we obtainnextand we easily see that the real part of this expression is negative.

Next, we calculateEventually,and hence,For our reference values the expression in the brackets equals .

At the end we sum up the results of the Hopf bifurcation analysis in the following corollary.

Corollary 5. For reference parameter values the system described by (1) undergoes two subsequent subcritical Hopf bifurcations.

4. Numerical Simulations

In order to simulate the decision-making process, we assume a transient change in one of the inputs to the nodes; that is, we start from the resting state of the network and then solve (1) with andwhere is the stimulation duration time. In a properly working decision-making network, we should expect that the bigger the value of is, that is, the stimulus, the faster the decision in favor of population according to the psychometric function (2) will be.

In all the numerical experiments, we have chosen and (for this value the time scale allows to reproduce delays present in real systems). Note that the qualitative dynamics of the system described by (1) does not depend on . Moreover, we choose a reference value to reflect that there is a certain baseline input to the network. Other parameter values fixed across all simulations were , , and  s. We perform numerical simulations of the model with respect to the magnitude of delay and stimulus strength . For all of the simulations we use the built-in MATLAB delayed differential equation solver dde23 [28] with lowered default tolerances (RelTol and AbsTol equal to ). We performed the simulations on the Neuroscience Gateway platform [29].

One important aspect of the considered (1) is that we cannot guarantee the nonnegativity of solutions as the delayed self-inhibition term has a negative sign. It is obvious, however, that the firing rate cannot have a negative value and, in order for the model to be physiologically valid, we need to impose a barrier for the firing rates at zero. Thus, in the simulations, whenever one of the coordinates reaches value of zero, we stop the simulation and solve the reduced system with one of the coordinates equal to zero until its derivative becomes positive again; that is, its delayed value becomes lower than .

We start our numerical experiments from simulations in which stimulus of a fixed magnitude is applied to the network with different delay in the self-inhibition terms. As it can be expected from the analysis, for small delays, that is, below the Hopf bifurcation threshold, the solutions exhibit vanishing oscillations around the stable steady state; see left panel in Figure 3(a). The stimulus in this particular case is too small for any certain decision to be made, and both psychometric functions remain separated from the value of 1; see left panel in Figure 3(b). It is clear, however, that if the stimulus was stronger, then a certain decision would be made, that is, would cross the threshold of at some point.

Figure 3: Comparison of the model solutions ((a); (1)) together with the corresponding psychometric function values ((b); (2) and (3)) for different values of time delay and stimulus strength . For small delays, that is, in the stability regime, we observe vanishing oscillations with no certain decision being made (). Interestingly, for which is above critical destabilizing delay value we observe multiple preference switches (middle column). For certain decision is being made ().

Interestingly, for delays above the critical value at which first Hopf bifurcation occurs, we observe very nonintuitive behavior; see middle panels in Figures 3(a) and 3(b). Namely, the psychometric function value shows multiple perceptual switches; that is, there are multiple time points in which there is a change from into . For even larger delays we observe that a certain decision has been made for the same stimulus of magnitude that was insufficient to make a certain decision in the case of small delays; see right panels Figures 3(a) and 3(b).

Because of the observed nonintuitive behavior of the psychometric function, we decided to calculate a decision map; that is, to evaluate the psychometric function values on the model solutions for different values of delay and the stimulus strength , see Figure 4. As expected, for delay values below the critical Hopf bifurcation threshold, the certainty of the decision depends directly on the strength of the stimulus, and the networks is always able to correctly identify that the stimulus was applied to population . Interestingly, this is not the case for larger delays and we observe that for the same value of delay above the critical first bifurcation threshold the decision can be opposite depending on the stimulus strength, compare Figure 4(a). Moreover, above the bifurcation threshold we can have more than 10 transient decision switches before any certain decision is made; compare Figure 4(b).