Applications of Delay Differential Equations in Biological SystemsView this Special Issue
Impact of Time Delay in Perceptual Decision-Making: Neuronal Population Modeling Approach
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.
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 . 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 [6–8]. In particular, the working memory—which involves active manipulation of information—is affected, and this effect can further influence the perceptual decision-making . (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 , 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  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 [14–16], 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 .
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  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 . 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  and computational  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).
In order to describe the temporal dynamics of considered system (Figure 1), we introduce a system of delay differential equations (DDEs ) 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 ; that is, the dynamics of synaptic weights is firing rate-dependent. For in-depth analytical and numerical investigations we have chosen a biologically plausible  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
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.
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).
Lemma 2 holds under weaker assumptions, that is, for and being analytic functions, not only polynomials; compare  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 .
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  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 . 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 . 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,