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

Authors: Urszula Foryś; Natalia Z. Bielczyk; Katarzyna Piskała; Martyna Płomecka; Jan Poleszczuk
Journal: Complexity
 (2017)

---

## Abstract

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.


---

## Body

## 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 [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 [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 [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 [14].The paper is organized in the following way. In Section2, 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 ratesr(t) of two positively interconnected self-inhibiting neuronal populations that receive external inputs I1(t) and I2(t), 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 I1(t) and I2(t), respectively. The dynamics of this system is described by (1).In order to describe the temporal dynamics of considered system (Figure1), we introduce a system of delay differential equations (DDEs [21]) in the following form:(1)r˙1t=αI1t-r1t-τ+ϵfr1tr2tr2t,r˙2t=αI2t-r2t-τ+ϵfr1tr2tr1t,where α=1/τr, τr 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 f 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 f(x)=x2/(1+x2).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 inputsI1(t)=I2(t)≡I 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 decodingwhich 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 r1(t)-r2(t), the more likely the decision that the population 1 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 of(2)p1t=11+exp⁡-β∫0tr1ξ-r2ξdξ,where β is a parameter influencing the slope of the sigmoid function with respect to the cumulative difference ∫0tr1(ξ)-r2(ξ)dξ. Similarly, we define the psychometric function for the second population as(3)p2t=11+exp⁡-β∫0tr2ξ-r1ξdξ.Values p1(t) and p2(t) can only asymptotically approach 1; therefore we add a condition that if, at a given time point t, p1(t) exceeds the threshold value of 1-γ (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 α=1 and constant symmetric input I. Notice that the qualitative dynamics of (1) does not depend on α.


### 3.1. Behavior of the Model forτ=0 and Constant Symmetric Input I

In this subsection, we present a detailed analysis of the model dynamics forτ=0 and I1(t)=I2(t)≡I, 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 equations(4)I-r1+ϵr2fr1r2=0,I-r2+ϵr1fr1r2=0,which yields r1-r2+ϵ(r1-r2)f(r1r2)=0. It is then obvious that both coordinates of any steady state are the same. Let us denote a steady state by (r¯,r¯). Clearly, I-r¯+ϵr¯f(r¯2)=0, and the number of steady states depends on the shape of the graph of hϵ(r)=(1-ϵ)r+ϵ(r/(1+r4)). Notice that the reference case ϵ=1 is specific, as for ϵ≠1 the function h is asymptotically linear ~(1-ϵ)r, while for ϵ=1 it tends to 0 as r→∞.Let us considerϵ=1. Then we have h1′(r)=(1-3r4)/1+r42, and therefore h1 is increasing for r∈[0,274/3], achieves its maximum h1m=274/4≈0.5699 at rm=274/3, and decreases to 0 for r>rm. This implies that there are two steady states for 0<I<h1m and no steady state for I>h1m, while for I=h1m there is a bifurcation. The reference value I=0.4<h1m, so there exist two steady states, r¯≈0.4115 and r¯¯≈1.1827.Forϵ≠1, we have hϵ′(r)=1-ϵ+ϵ((1-3r4)/1+r42), and looking for zeros of hϵ′ we obtain (1-ϵ)r8+(2-5ϵ)r4+1=0, and we see that this quadratic equation has two positive solutions for ϵ∈(16/25,1), no real solution for ϵ<16/25, and one positive solution for ϵ>1. This means that for ϵ<16/25 the function hϵ is increasing, for ϵ≥1 it has one maximum and tends either to 0 (for ϵ=1) or to -∞ (for ϵ>1), while for ϵ∈(16/25,1) it is first increasing, then decreasing, and eventually increasing linearly to +∞.Corollary 1.
Forϵ<16/25 there is one steady state of (1); for ϵ≥1 there are two steady states for small values of I and no steady state for larger I values; for ϵ∈(16/25,1) there is one steady state for small and sufficiently large values of I, while for intermediate I 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 r2=r1. Clearly, assuming r1=r2=r from both equations of (1) we obtain(5)r˙=I-r+ϵrfr2.The number of steady states determines the dynamics of (5). Let us assume that ϵ<16/25. Then there is only one steady state r¯ and the right-hand side of (5) is positive for r<r¯ and negative for r>r¯. Therefore, r¯ is globally attractive. If ϵ≥1, then there are two steady states r¯<r¯¯ for small I, and the right-hand side is positive for r<r¯ and r>r¯¯ and negative for r∈(r¯,r¯¯). This means that r¯ is locally stable, while r¯¯ is unstable. Moreover, solutions for r(0)>r¯¯ tends to ∞. If there is no steady state, then all solutions tend to +∞, as they are increasing and unbounded. For ϵ∈(16/25,1) there can be up to three steady states. Assume that there are three steady states r¯<r¯¯<r¯¯¯. In this case r¯ and r¯¯¯ are stable, while r¯¯ is unstable. Moreover, for r(0)>r¯¯¯ solutions tend to r¯¯¯.It is obvious that the symmetry of (1) implies that the phase space is divided into two symmetric subspaces by the straight line r2=r1. Moreover, the dynamics of (5) is crucial in determining the whole model dynamics. Notice that r˙1=I-r1+ϵr2f(r1r2)<I-r1+ϵr2 and for r2<(r1-I)/ϵ we have r˙1<0. Similarly, for r2>I+ϵr1 we have r˙2<0. Moreover, these straight lines are asymptotes for null-clines of (1). Therefore, asymptotic dynamics for large values of r1, r2 is related to the dynamics of (5) for large r, 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 I=0.4 and ϵ=0.5 (a) (there is only one steady state, stable node at around (0.41,0.41)); ϵ=0.87 (b) (there are three steady states: stable nodes at around (0.41,0.41) and (2.77,2.77) and a saddle at around (1.55,1.55)); ϵ=1 (c) (there are two steady states: a stable node at around (0.41,0.41) and a saddle at around (1.18,1.18)).
(a)
(b)
(c)


### 3.2. Model Behavior forτ>0 and Constant Symmetric Input I

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τth>0 for which the steady state (r¯,r¯) 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 DDE(6)x˙=gxt,g0=0,with a smooth function g. Let W(λ,τ) denote a characteristic function for (6) at x¯=0. Assume that(7)Wλ,τ=WIλ,τ·WIIλ,τ,WIλ,τ=Pλ+Qλe-λτ,where P and Q are polynomials, deg⁡P>deg⁡Q. Together with (6) we consider(8)x˙=g1xt,g10=0,for which WI is a characteristic function.Lemma 2.
Assume thatP and Q have no common imaginary root and WI has a pair of purely imaginary simple eigenvalues ±iω0(ω0>0) for some critical value τ=τ0>0, and moreover these eigenvalues satisfy transversality condition d/dτRe⁡λ(τ)τ=τ0≠0. If WII(iω0,τ)≠0, then ±iω0 is a pair of simple eigenvalues of W(λ,τ) for τ=τ0 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(d/dτ)Re⁡λ(τ)τ=τ0 for (6) is obtained from the relation(9)ddτWIλ,ττ=τ0·WIIλ,ττ=τ0+WIλ,ττ=τ0·ddτWIIλ,ττ=τ0=0.As WI(λ,τ)|τ=τ0=0 while WII(λ,τ)|τ=τ0≠0, this implies (d/dτ)(WI(λ,τ))τ=τ0≠0 which is the relation determining transversality condition for (8).Lemma2 holds under weaker assumptions, that is, for P and Q 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.
LetWI(λ,τ)=λ+γ+e-λτ, |γ|<1. Then there exists a pair of purely imaginary eigenvalues ±iω0, ω0=1-γ2, for τ0=arccos⁡(-γ)/ω0  for which eigenvalues cross imaginary axis from left to right. Moreover, a steady state x¯=0 loses stability for τ=τ0 and cannot gain it again for τ>τ0.Proof.
Following [24] we define an auxiliary function FA(ω)=|P(iω)|2-|Q(iω)|2=ω2+γ2-1. It is obvious that ω0=1-γ2 is a simple zero of FA. Looking for critical delay related to this pair of eigenvalues ±iω0, we need to solve a system of equations(10)cos⁡ω0τ=-γ,sin⁡ω0τ=ω0>0.It is obvious that we obtain a sequence of critical delays τn=(1/ω0)arccos⁡(-γ)+2nπ, n∈N. However, as the derivative of the auxiliary function is positive, eigenvalues always cross imaginary axis from left to right (independently of n), which means that a switch of stability appears only at τ=τ0.Now, we turn to the main topic of this subsection that is analysis of (1) for τ>0. 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,(11)X˙=LXt+GXt,where X=(x,y)=(r1-r¯,r2-r¯) and Xt(h)=X(t+h) for h∈[-τ,0], L is a linear part, and G is nonlinear part of the system. Hence,(12)Lxt,yt=L1xt,yt,L2xt,yt=-xt-τ+ηxt0+η+βyt0,-yt-τ+ηyt0+η+βxt0T,Gxt,yt=G1xt,yt,G2xt,yt=yt0+r¯Fxt,yt-β-ηxt0+yt0,xt0+r¯Fxt,yt-β-ηxt0+yt0T,where F(xt,yt)=F(x(t),y(t))=ϵf((x(t)+r¯)(y(t)+r¯)), η=ϵf′(r¯2)r¯2>0, and β=ϵf(r¯2)>0. With the linear operator L we are able to associate an operator T(t)X0=Xt which is a solution to (1) for initial data X0. In this way we obtain a strongly continuous semigroup generated by an infinitesimal generator A (cf. [26, 27]). If for some critical value τ0 the generator A has an eigenvalue λ=iω0, then a Hopf bifurcation occurs under the assumptions that the eigenvalues ±iω0 are simple and cross imaginary axis with nonzero speed when τ crosses τ0.Let us denote a characteristic matrix byΔ(λ,τ); that is,(13)Δλ,τ=λ+exp⁡-λτ-η-β+η-β+ηλ+exp⁡-λτ-η.Looking for eigenvalues, we need to find zeros of the characteristic function(14)Wλ,τ=det⁡Δλ,τ=λ+exp⁡-λτ-η2-η+β2.It is obvious that(15)Wλ,τ=WIλ,τ·WIIλ,τ,WIλ,τ=λ-2η-β+e-λτ,WIIλ,τ=λ+β+e-λτ.Clearly, we can use Lemmas 2 and 3 in the analysis of stability switches for the steady state (r¯,r¯). In the considered case, the steady state (r¯,r¯) is stable for τ=0; that is, W(λ,0)=λ-2η-β+1λ+β+1 has negative zeros, yielding 1>2η+β. Hence, both quasi-polynomials WI and WII satisfy assumptions of Lemma 3. This means that there are two sequences of critical delays τnIn∈N and τnIIn∈N associated with WI and WII, respectively. However, the switch of stability can occur only for τ0I or τ0II, depending on the magnitude of those delays. Clearly, according to Lemma 2 eigenvalues cross imaginary axis in the same direction for both WI and WII, which means that they cross from left to right, and therefore the steady state (r¯,r¯) loses stability for the smallest critical delay.As a result, we can state the following theorem.Theorem 4.
Letτ0=arccos⁡(2η+β)/1-(2η+β)2. The steady state (r¯,r¯) of (1) is locally asymptotically stable for τ<τ0 and unstable for τ>τ0, and at τ=τ0 a Hopf bifurcation occurs.Proof.
We only need to check thatτ0=τ0I<τ0II. Notice that τ0II=arccos⁡(-β)/1-β2. Moreover, the function arccos⁡(x)/1-x2, x∈(-1,1) is decreasing. As arccos⁡(-β)>arccos⁡(β) we obtain(16)τ0II=arccos⁡-β1-β2>arccos⁡β1-β2>arccos⁡2η+β1-2η+β2=τ0I.For the reference values of parameters forr¯≈0.4115 we obtain β≈0.02787 and η≈0.05418. Considering WI we obtain ω0I≈0.9907 and τ0I≈1.4476. For WII we have ω0II≈0.9996 and τ0II=arccos⁡(-β)/1-β2≈1.5993.Hence, we study the bifurcation atτ0=τ0I and corresponding ω0=1-(2η+β)2. Let us denote β1=2η+β to shorten the notation. Therefore, ω0=1-β12. Notice that from the form of WI we obtain iω0+e-iω0τ0=β1.In the following, we base on the ideas presented in [25]. We know that λ=iω0 is a purely imaginary eigenvalue of the infinitesimal generator A for τ=τ0 if there exists a vector p∈C2 such that Δ(iω0,τ0)p=0 and then Φ(h)=eiω0hp is an eigenvector for A at τ0. Moreover, iω0 is also an eigenvalue for the adjoint operator A∗ for τ=τ0 and Ψ(h)=qeiω0h is an eigenvector, where qΔ(iω0,τ0)=0. In the considered case, we are able to choose q such that qd1Δ(iω0,τ0)p=1, where d1 is the derivative with respect to the first variable, here λ.When looking forp=(a,b)T we obtain(17)Δλ,τp=λ+exp⁡-λτ-η-β+η-β+ηλ+exp⁡-λτ-η·ab=η+β-η+β-η+βη+β·ab=0,0T,and therefore(18)β+ηa-β+ηb=0.Clearly, a=b and we can choose p=1,1T. Hence, Φ(h)=eiω0h(1,1)T is the eigenvector for the eigenvalue iω0.Next, we need to find a vectorq that satisfies(19)qΔλ,τ=a,b·η+β-η+β-η+βη+β=0,0,which means that coordinates of q satisfy the same relation as for p, that is, a=b. Next, calculating(20)d1Δiω0,τ0=1-τ0exp⁡-iω0τ0001-τ0exp⁡-iω0τ0we obtain 2a(1-τ0exp⁡(-iω0τ0))=1, that is, a=1/21+τ0(iω0-β1), and eventually(21)a=121-τ0β1-iω0τ01-τ0β12+ω02τ02.The type of the studied bifurcation is determined by a coefficientμ2 of the third term in Taylor expansion of an orbit. This coefficient reads(22)μ2=Re⁡cRe⁡qd2Δiω0,τ0p,where d2 denotes a first derivative with respect to the second variable τ, while c=cI+cII+cIII, and(23)cI=12qd13G0,τ0Φ,Φ,Φ¯,cII=qd12G0,τ0Ψ1,Φ,cIII=12qd12G0,τ0Ψ2,Φ¯,where d1i,i=2,3, denotes ith derivative with respect to the first variable, Ψ1(h)=Δ(0,τ0)-1d12G(0,τ0)Φ,Φ¯ (in fact Ψ1 is a constant function as it does not depend on h), and Ψ2(h)=e2iω0hΔ(2iω0,τ0)-1d12G(0,τ0)Φ,Φ.Ifμ2>0, then the bifurcation is supercritical; that is, periodic solutions appear for τ>τ0 and are stable in such a case. If μ2<0, then the bifurcation is subcritical: periodic solutions exist for τ<τ0, and as the steady state is stable in this case, periodic orbits are necessarily unstable.First, we calculate the denominator ofμ2. We have(24)d2Δiω0,τ0=-iω0exp⁡-iω0τ000-iω0exp⁡-iω0τ0=-ω02-iω0β100-ω02-iω0β1,q·d2Δiω0,τ0·p=-2aω02+iω0β1,and we easily check that the real part of this expression is equal to -ω02/(1-τ0β12+ω02τ02)<0.Next, we would like to calculate the numerator ofμ2. To this end, we need to calculate the derivatives of the nonlinear part G necessary to calculate c, and we omit bar in the notation r¯. Denoting by u, v, and w test functions from C[-τ,0],R2, we obtain(25)d1G1xt,ytu=yt+rFx′xt,yt-ηu10+Fxt,yt+yt+rFy′xt,yt-β-ηu20,d1G2xt,ytu=Fxt,yt+xt+rFx′xt,yt-β-ηu10+xt+rFy′xt,yt-ηu20,d12G1xt,ytu,v=yt+rFxx′′xt,ytu10v10+Fx′xt,yt+yt+rFxy′′xt,ytu10v20+u20v10+2Fy′xt,yt+yt+rFyy′′xt,ytu20v20,d12G2xt,ytu,v=2Fx′xt,yt+xt+rFxx′′xt,ytu10v10+Fy′xt,yt+xt+rFxy′′xt,ytu10v20+u20v10+xt+r¯Fyy′′xt,ytu20v20,d13G1xt,ytu,v,w=yt+rFxxx′′′xt,ytu10v10w10+Fxx′′xt,yt+yt+rFxxy′′′xt,ytu10v10w20+u10v20w10+u20v10w10+2Fxy′′xt,yt+yt+rFxyy′′′xt,ytu10v20w20+u20v10w20+u20v20w10+3Fyy′′xt,yt+yt+rFyyy′′′xt,ytu20v20w20,d13G2xt,ytu,v,w=3Fxx′′xt,yt+xt+rFxxx′′′xt,ytu10v10w10+2Fxy′′xt,yt+xt+rFxxy′′′xt,ytu10v10w20+u10v20w10+u20v10w10+Fyy′′xt,yt+xt+rFxyy′′′xt,ytu10v20w20+u20v10w20+u20v20w10+xt+rFyyy′′′xt,ytu20v20w20.Evaluating the second and third derivative at (0,0) we obtain(26)d12G1u,v=2ϵr31+r431-3r4u10v10+3-r4u10v20+u20v10+3-r4u20v20,d12G2u,v=2ϵr31+r433-r4u10v10+3-r4u10v20+u20v10+1-3r4u20v20,d13G1u,v,w=2ϵr21+r4412r4r4-1u10v10w10+31-6r4+r8u10v10w20+u10v20w10+u20v10w10+23-8r4+r8u10v20w20+u20v10w20+u20v20w10+31-6r4+r8u20v20w20,d13G2u,v,w=2ϵr21+r4431-6r4+r8u10v10w10+23-8r4+r8u10v10w20+u10v20w10+u20v10w10+31-6r4+r8u10v20w20+u20v10w20+u20v20w10+12r4r4-1u20v20w20,where d1jGi(0,0),i=1,2,j=2,3, is denoted by d1jGi, to shorten the notation.Using the formula above, we calculate(27)cI=12q·d13GΦ,Φ,Φ¯=12a,a·12ϵr25r8-22r4+51+r4412ϵr25r8-22r4+51+r44=12ϵr25r8-22r4+51+r44a.To calculate cII we need to evaluate (Δ(0,τ0))-1. We have(28)Δ0,τ0-1=11+β1-β11-ηη+βη+β1-η.Moreover,(29)d12GΦ,Φ¯=4ϵr35-3r41+r434ϵr35-3r41+r43,and therefore(30)Ψ1=Δ0,τ0-1d12GΦ,Φ¯=4ϵr35-3r41+r431-β111.Finally,(31)cII=q·d12GΨ1,Φ=4ϵr35-3r41+r431-β1u,u·4ϵr35-3r41+r434ϵr35-3r41+r43=32ϵ2r65-3r421+r461-β1u.To obtain the last parameter cIII we first need to calculate (Δ(2iω0,τ0))-1.Let us denoteM=Δ(2iω0,τ0) and calculate the first term of this matrix; that is,(32)2iω0+e-2iω0τ0-η=2iω0+β1-iω02-η=β12-ω02-η+2iω01-β1.Hence,(33)det⁡M=β12-ω02-η+2iω01-β12-η+β2=β12-ω02+β+2iω01-β1β12-ω02-β1+2iω01-β1,M-1=1det⁡Mβ12-ω02-η+2iω01-β1η+βη+ββ12-ω02-η+2iω01-β1.Next we calculate Ψ2(h)=e2iω0hV, where(34)V=M-1d12GΦ,Φ=M-14ϵr35-3r41+r434ϵr35-3r41+r43=4ϵr35-3r4/1+r43β12-ω02+β+2iω01-β1det⁡M11=4ϵr35-3r41+r43β12-ω02-β1+2iω01-β111.Finally, we evaluate(35)d12GΨ2,Φ¯=4ϵr35-3r41+r43β12-ω02-β1+2iω01-β14ϵr35-3r41+r434ϵr35-3r41+r43=16ϵ2r65-3r421+r46β12-ω02-β1+2iω01-β111,and therefore(36)cIII=12a,a·d12GΨ2,Φ¯=16ϵ2r65-3r421+r46β12-ω02-β1+2iω01-β1a.Now, we are in a position to check the sign ofRe⁡c, where c=cI+cII+cIII. Clearly,(37)c=4ϵr21+r4435r8-22r4+5+8ϵr45-3r421+r421-β1+4ϵr45-3r42β12-ω02-β1-2iω01-β11+r42β12-ω02-β12+4ω021-β12a,and therefore(38)1+r444ϵr2Re⁡c=35r8-22r4+5+8ϵr45-3r421+r421-β1+4ϵr45-3r42β12-ω02-β11+r42β12-ω02-β12+4ω021-β12Re⁡a+8ω0ϵr45-3r421-β11+r42β12-ω02-β12+4ω021-β12Im⁡a,where Re⁡a=(1-τ0β1)/21-τ0β12+ω02τ02 and Im⁡a=-ω0τ0/21-τ0β12+ω02τ02.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 r¯<0.5, the first term is positive and of the order of units, suggesting that Re⁡c>0. 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 Re⁡c, as it is difficult to guess its sign.Now, we calculateRe⁡c for the reference parameters with ϵ=1. Let us recall that r=0.4115, β=0.02787, η=0.05418, ω0=0.9907, τ0=1.4476, and moreover β1=2η+β=0.13623. We obtain Re⁡a≈0.1486 and Im⁡a≈-0.2655. Consecutive terms in the brackets are equal to 13.11988735, 6.060143781, and -0.6953365688, while the last fraction is 1.082692850. Eventually,(39)Re⁡c≈0.604908396418.48469456·0.1486-1.082692850·0.2655≈1.487693962>0,which means that, in our reference case, the bifurcation issubcritical.In order to complement the analysis above, we can also show that in the reference case the second bifurcation appearing forτ0II is subcritical as well. We again denote τ0II by τ0, to shorten the notation. Now, the pair ω0, τ0 satisfies iω0+β+exp⁡(-iω0τ0)=0 yielding iω0+exp⁡(-iω0τ0)=-β, and from this relation we obtain the vector p=(1,-1)T. Hence, Φ(h)=exp⁡(iω0h)(1,-1)T is the eigenvector for the eigenvalue iω0. Next, we find q such that qd1Δ(iω0,τ0)p=1. As before, coordinates of q satisfies the same relation as for p, that is, q=(a,a), a==(1/2)((1+τ0β-iω0τ0)/((1+τ0β)2+ω02τ02)). Calculating the denominator of μ2 we obtain(40)d2Δiω0,τ0=-ω02+iω0β00-ω02+iω0β,next(41)q·d2Δiω0,τ0·p=2a-ω02+iω0β=iω0τ0-τ0β+1ω02τ02+τ0β+12ω02-iω0β=-ω0+iβ+τ0β2+ω02ω02τ02+τ0β+12ω0,and we easily see that the real part of this expression is negative.Next, we calculate(42)cI=12q·d13GΦ,Φ,Φ¯=12a,-a·12ϵr21+r42-12ϵr21+r42=12ϵr21+r42a,cII=-4ϵr31+r4211-β1a,-a·-4ϵr31+r424ϵr31+r42=4ϵr321+r442a1-β1,cIII=4r3ϵ2u1+r442iω0β+1+β2-ω02-β1.Eventually,(43)c=4ϵr21+r4431+r42+8ϵr41-β1+4ϵr4β2-ω02-β1-2iω0β+14ω02β+12+β2-ω02β12a,and hence,(44)sign⁡Re⁡c=sign⁡1+τ0β31+r42+8ϵr41-β1+4ϵr4β2-ω02-β14ω02β+12+β2-ω02-β12-8ϵr2τ0ω02β+14ω02β+12+β2-ω02-β12.For our reference values the expression in the brackets equals 3.146913951>0.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.


## 3.1. Behavior of the Model forτ=0 and Constant Symmetric Input I

In this subsection, we present a detailed analysis of the model dynamics forτ=0 and I1(t)=I2(t)≡I, 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 equations(4)I-r1+ϵr2fr1r2=0,I-r2+ϵr1fr1r2=0,which yields r1-r2+ϵ(r1-r2)f(r1r2)=0. It is then obvious that both coordinates of any steady state are the same. Let us denote a steady state by (r¯,r¯). Clearly, I-r¯+ϵr¯f(r¯2)=0, and the number of steady states depends on the shape of the graph of hϵ(r)=(1-ϵ)r+ϵ(r/(1+r4)). Notice that the reference case ϵ=1 is specific, as for ϵ≠1 the function h is asymptotically linear ~(1-ϵ)r, while for ϵ=1 it tends to 0 as r→∞.Let us considerϵ=1. Then we have h1′(r)=(1-3r4)/1+r42, and therefore h1 is increasing for r∈[0,274/3], achieves its maximum h1m=274/4≈0.5699 at rm=274/3, and decreases to 0 for r>rm. This implies that there are two steady states for 0<I<h1m and no steady state for I>h1m, while for I=h1m there is a bifurcation. The reference value I=0.4<h1m, so there exist two steady states, r¯≈0.4115 and r¯¯≈1.1827.Forϵ≠1, we have hϵ′(r)=1-ϵ+ϵ((1-3r4)/1+r42), and looking for zeros of hϵ′ we obtain (1-ϵ)r8+(2-5ϵ)r4+1=0, and we see that this quadratic equation has two positive solutions for ϵ∈(16/25,1), no real solution for ϵ<16/25, and one positive solution for ϵ>1. This means that for ϵ<16/25 the function hϵ is increasing, for ϵ≥1 it has one maximum and tends either to 0 (for ϵ=1) or to -∞ (for ϵ>1), while for ϵ∈(16/25,1) it is first increasing, then decreasing, and eventually increasing linearly to +∞.Corollary 1.
Forϵ<16/25 there is one steady state of (1); for ϵ≥1 there are two steady states for small values of I and no steady state for larger I values; for ϵ∈(16/25,1) there is one steady state for small and sufficiently large values of I, while for intermediate I 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 r2=r1. Clearly, assuming r1=r2=r from both equations of (1) we obtain(5)r˙=I-r+ϵrfr2.The number of steady states determines the dynamics of (5). Let us assume that ϵ<16/25. Then there is only one steady state r¯ and the right-hand side of (5) is positive for r<r¯ and negative for r>r¯. Therefore, r¯ is globally attractive. If ϵ≥1, then there are two steady states r¯<r¯¯ for small I, and the right-hand side is positive for r<r¯ and r>r¯¯ and negative for r∈(r¯,r¯¯). This means that r¯ is locally stable, while r¯¯ is unstable. Moreover, solutions for r(0)>r¯¯ tends to ∞. If there is no steady state, then all solutions tend to +∞, as they are increasing and unbounded. For ϵ∈(16/25,1) there can be up to three steady states. Assume that there are three steady states r¯<r¯¯<r¯¯¯. In this case r¯ and r¯¯¯ are stable, while r¯¯ is unstable. Moreover, for r(0)>r¯¯¯ solutions tend to r¯¯¯.It is obvious that the symmetry of (1) implies that the phase space is divided into two symmetric subspaces by the straight line r2=r1. Moreover, the dynamics of (5) is crucial in determining the whole model dynamics. Notice that r˙1=I-r1+ϵr2f(r1r2)<I-r1+ϵr2 and for r2<(r1-I)/ϵ we have r˙1<0. Similarly, for r2>I+ϵr1 we have r˙2<0. Moreover, these straight lines are asymptotes for null-clines of (1). Therefore, asymptotic dynamics for large values of r1, r2 is related to the dynamics of (5) for large r, 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 I=0.4 and ϵ=0.5 (a) (there is only one steady state, stable node at around (0.41,0.41)); ϵ=0.87 (b) (there are three steady states: stable nodes at around (0.41,0.41) and (2.77,2.77) and a saddle at around (1.55,1.55)); ϵ=1 (c) (there are two steady states: a stable node at around (0.41,0.41) and a saddle at around (1.18,1.18)).
(a)
(b)
(c)


## 3.2. Model Behavior forτ>0 and Constant Symmetric Input I

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τth>0 for which the steady state (r¯,r¯) 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 DDE(6)x˙=gxt,g0=0,with a smooth function g. Let W(λ,τ) denote a characteristic function for (6) at x¯=0. Assume that(7)Wλ,τ=WIλ,τ·WIIλ,τ,WIλ,τ=Pλ+Qλe-λτ,where P and Q are polynomials, deg⁡P>deg⁡Q. Together with (6) we consider(8)x˙=g1xt,g10=0,for which WI is a characteristic function.Lemma 2.
Assume thatP and Q have no common imaginary root and WI has a pair of purely imaginary simple eigenvalues ±iω0(ω0>0) for some critical value τ=τ0>0, and moreover these eigenvalues satisfy transversality condition d/dτRe⁡λ(τ)τ=τ0≠0. If WII(iω0,τ)≠0, then ±iω0 is a pair of simple eigenvalues of W(λ,τ) for τ=τ0 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(d/dτ)Re⁡λ(τ)τ=τ0 for (6) is obtained from the relation(9)ddτWIλ,ττ=τ0·WIIλ,ττ=τ0+WIλ,ττ=τ0·ddτWIIλ,ττ=τ0=0.As WI(λ,τ)|τ=τ0=0 while WII(λ,τ)|τ=τ0≠0, this implies (d/dτ)(WI(λ,τ))τ=τ0≠0 which is the relation determining transversality condition for (8).Lemma2 holds under weaker assumptions, that is, for P and Q 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.
LetWI(λ,τ)=λ+γ+e-λτ, |γ|<1. Then there exists a pair of purely imaginary eigenvalues ±iω0, ω0=1-γ2, for τ0=arccos⁡(-γ)/ω0  for which eigenvalues cross imaginary axis from left to right. Moreover, a steady state x¯=0 loses stability for τ=τ0 and cannot gain it again for τ>τ0.Proof.
Following [24] we define an auxiliary function FA(ω)=|P(iω)|2-|Q(iω)|2=ω2+γ2-1. It is obvious that ω0=1-γ2 is a simple zero of FA. Looking for critical delay related to this pair of eigenvalues ±iω0, we need to solve a system of equations(10)cos⁡ω0τ=-γ,sin⁡ω0τ=ω0>0.It is obvious that we obtain a sequence of critical delays τn=(1/ω0)arccos⁡(-γ)+2nπ, n∈N. However, as the derivative of the auxiliary function is positive, eigenvalues always cross imaginary axis from left to right (independently of n), which means that a switch of stability appears only at τ=τ0.Now, we turn to the main topic of this subsection that is analysis of (1) for τ>0. 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,(11)X˙=LXt+GXt,where X=(x,y)=(r1-r¯,r2-r¯) and Xt(h)=X(t+h) for h∈[-τ,0], L is a linear part, and G is nonlinear part of the system. Hence,(12)Lxt,yt=L1xt,yt,L2xt,yt=-xt-τ+ηxt0+η+βyt0,-yt-τ+ηyt0+η+βxt0T,Gxt,yt=G1xt,yt,G2xt,yt=yt0+r¯Fxt,yt-β-ηxt0+yt0,xt0+r¯Fxt,yt-β-ηxt0+yt0T,where F(xt,yt)=F(x(t),y(t))=ϵf((x(t)+r¯)(y(t)+r¯)), η=ϵf′(r¯2)r¯2>0, and β=ϵf(r¯2)>0. With the linear operator L we are able to associate an operator T(t)X0=Xt which is a solution to (1) for initial data X0. In this way we obtain a strongly continuous semigroup generated by an infinitesimal generator A (cf. [26, 27]). If for some critical value τ0 the generator A has an eigenvalue λ=iω0, then a Hopf bifurcation occurs under the assumptions that the eigenvalues ±iω0 are simple and cross imaginary axis with nonzero speed when τ crosses τ0.Let us denote a characteristic matrix byΔ(λ,τ); that is,(13)Δλ,τ=λ+exp⁡-λτ-η-β+η-β+ηλ+exp⁡-λτ-η.Looking for eigenvalues, we need to find zeros of the characteristic function(14)Wλ,τ=det⁡Δλ,τ=λ+exp⁡-λτ-η2-η+β2.It is obvious that(15)Wλ,τ=WIλ,τ·WIIλ,τ,WIλ,τ=λ-2η-β+e-λτ,WIIλ,τ=λ+β+e-λτ.Clearly, we can use Lemmas 2 and 3 in the analysis of stability switches for the steady state (r¯,r¯). In the considered case, the steady state (r¯,r¯) is stable for τ=0; that is, W(λ,0)=λ-2η-β+1λ+β+1 has negative zeros, yielding 1>2η+β. Hence, both quasi-polynomials WI and WII satisfy assumptions of Lemma 3. This means that there are two sequences of critical delays τnIn∈N and τnIIn∈N associated with WI and WII, respectively. However, the switch of stability can occur only for τ0I or τ0II, depending on the magnitude of those delays. Clearly, according to Lemma 2 eigenvalues cross imaginary axis in the same direction for both WI and WII, which means that they cross from left to right, and therefore the steady state (r¯,r¯) loses stability for the smallest critical delay.As a result, we can state the following theorem.Theorem 4.
Letτ0=arccos⁡(2η+β)/1-(2η+β)2. The steady state (r¯,r¯) of (1) is locally asymptotically stable for τ<τ0 and unstable for τ>τ0, and at τ=τ0 a Hopf bifurcation occurs.Proof.
We only need to check thatτ0=τ0I<τ0II. Notice that τ0II=arccos⁡(-β)/1-β2. Moreover, the function arccos⁡(x)/1-x2, x∈(-1,1) is decreasing. As arccos⁡(-β)>arccos⁡(β) we obtain(16)τ0II=arccos⁡-β1-β2>arccos⁡β1-β2>arccos⁡2η+β1-2η+β2=τ0I.For the reference values of parameters forr¯≈0.4115 we obtain β≈0.02787 and η≈0.05418. Considering WI we obtain ω0I≈0.9907 and τ0I≈1.4476. For WII we have ω0II≈0.9996 and τ0II=arccos⁡(-β)/1-β2≈1.5993.Hence, we study the bifurcation atτ0=τ0I and corresponding ω0=1-(2η+β)2. Let us denote β1=2η+β to shorten the notation. Therefore, ω0=1-β12. Notice that from the form of WI we obtain iω0+e-iω0τ0=β1.In the following, we base on the ideas presented in [25]. We know that λ=iω0 is a purely imaginary eigenvalue of the infinitesimal generator A for τ=τ0 if there exists a vector p∈C2 such that Δ(iω0,τ0)p=0 and then Φ(h)=eiω0hp is an eigenvector for A at τ0. Moreover, iω0 is also an eigenvalue for the adjoint operator A∗ for τ=τ0 and Ψ(h)=qeiω0h is an eigenvector, where qΔ(iω0,τ0)=0. In the considered case, we are able to choose q such that qd1Δ(iω0,τ0)p=1, where d1 is the derivative with respect to the first variable, here λ.When looking forp=(a,b)T we obtain(17)Δλ,τp=λ+exp⁡-λτ-η-β+η-β+ηλ+exp⁡-λτ-η·ab=η+β-η+β-η+βη+β·ab=0,0T,and therefore(18)β+ηa-β+ηb=0.Clearly, a=b and we can choose p=1,1T. Hence, Φ(h)=eiω0h(1,1)T is the eigenvector for the eigenvalue iω0.Next, we need to find a vectorq that satisfies(19)qΔλ,τ=a,b·η+β-η+β-η+βη+β=0,0,which means that coordinates of q satisfy the same relation as for p, that is, a=b. Next, calculating(20)d1Δiω0,τ0=1-τ0exp⁡-iω0τ0001-τ0exp⁡-iω0τ0we obtain 2a(1-τ0exp⁡(-iω0τ0))=1, that is, a=1/21+τ0(iω0-β1), and eventually(21)a=121-τ0β1-iω0τ01-τ0β12+ω02τ02.The type of the studied bifurcation is determined by a coefficientμ2 of the third term in Taylor expansion of an orbit. This coefficient reads(22)μ2=Re⁡cRe⁡qd2Δiω0,τ0p,where d2 denotes a first derivative with respect to the second variable τ, while c=cI+cII+cIII, and(23)cI=12qd13G0,τ0Φ,Φ,Φ¯,cII=qd12G0,τ0Ψ1,Φ,cIII=12qd12G0,τ0Ψ2,Φ¯,where d1i,i=2,3, denotes ith derivative with respect to the first variable, Ψ1(h)=Δ(0,τ0)-1d12G(0,τ0)Φ,Φ¯ (in fact Ψ1 is a constant function as it does not depend on h), and Ψ2(h)=e2iω0hΔ(2iω0,τ0)-1d12G(0,τ0)Φ,Φ.Ifμ2>0, then the bifurcation is supercritical; that is, periodic solutions appear for τ>τ0 and are stable in such a case. If μ2<0, then the bifurcation is subcritical: periodic solutions exist for τ<τ0, and as the steady state is stable in this case, periodic orbits are necessarily unstable.First, we calculate the denominator ofμ2. We have(24)d2Δiω0,τ0=-iω0exp⁡-iω0τ000-iω0exp⁡-iω0τ0=-ω02-iω0β100-ω02-iω0β1,q·d2Δiω0,τ0·p=-2aω02+iω0β1,and we easily check that the real part of this expression is equal to -ω02/(1-τ0β12+ω02τ02)<0.Next, we would like to calculate the numerator ofμ2. To this end, we need to calculate the derivatives of the nonlinear part G necessary to calculate c, and we omit bar in the notation r¯. Denoting by u, v, and w test functions from C[-τ,0],R2, we obtain(25)d1G1xt,ytu=yt+rFx′xt,yt-ηu10+Fxt,yt+yt+rFy′xt,yt-β-ηu20,d1G2xt,ytu=Fxt,yt+xt+rFx′xt,yt-β-ηu10+xt+rFy′xt,yt-ηu20,d12G1xt,ytu,v=yt+rFxx′′xt,ytu10v10+Fx′xt,yt+yt+rFxy′′xt,ytu10v20+u20v10+2Fy′xt,yt+yt+rFyy′′xt,ytu20v20,d12G2xt,ytu,v=2Fx′xt,yt+xt+rFxx′′xt,ytu10v10+Fy′xt,yt+xt+rFxy′′xt,ytu10v20+u20v10+xt+r¯Fyy′′xt,ytu20v20,d13G1xt,ytu,v,w=yt+rFxxx′′′xt,ytu10v10w10+Fxx′′xt,yt+yt+rFxxy′′′xt,ytu10v10w20+u10v20w10+u20v10w10+2Fxy′′xt,yt+yt+rFxyy′′′xt,ytu10v20w20+u20v10w20+u20v20w10+3Fyy′′xt,yt+yt+rFyyy′′′xt,ytu20v20w20,d13G2xt,ytu,v,w=3Fxx′′xt,yt+xt+rFxxx′′′xt,ytu10v10w10+2Fxy′′xt,yt+xt+rFxxy′′′xt,ytu10v10w20+u10v20w10+u20v10w10+Fyy′′xt,yt+xt+rFxyy′′′xt,ytu10v20w20+u20v10w20+u20v20w10+xt+rFyyy′′′xt,ytu20v20w20.Evaluating the second and third derivative at (0,0) we obtain(26)d12G1u,v=2ϵr31+r431-3r4u10v10+3-r4u10v20+u20v10+3-r4u20v20,d12G2u,v=2ϵr31+r433-r4u10v10+3-r4u10v20+u20v10+1-3r4u20v20,d13G1u,v,w=2ϵr21+r4412r4r4-1u10v10w10+31-6r4+r8u10v10w20+u10v20w10+u20v10w10+23-8r4+r8u10v20w20+u20v10w20+u20v20w10+31-6r4+r8u20v20w20,d13G2u,v,w=2ϵr21+r4431-6r4+r8u10v10w10+23-8r4+r8u10v10w20+u10v20w10+u20v10w10+31-6r4+r8u10v20w20+u20v10w20+u20v20w10+12r4r4-1u20v20w20,where d1jGi(0,0),i=1,2,j=2,3, is denoted by d1jGi, to shorten the notation.Using the formula above, we calculate(27)cI=12q·d13GΦ,Φ,Φ¯=12a,a·12ϵr25r8-22r4+51+r4412ϵr25r8-22r4+51+r44=12ϵr25r8-22r4+51+r44a.To calculate cII we need to evaluate (Δ(0,τ0))-1. We have(28)Δ0,τ0-1=11+β1-β11-ηη+βη+β1-η.Moreover,(29)d12GΦ,Φ¯=4ϵr35-3r41+r434ϵr35-3r41+r43,and therefore(30)Ψ1=Δ0,τ0-1d12GΦ,Φ¯=4ϵr35-3r41+r431-β111.Finally,(31)cII=q·d12GΨ1,Φ=4ϵr35-3r41+r431-β1u,u·4ϵr35-3r41+r434ϵr35-3r41+r43=32ϵ2r65-3r421+r461-β1u.To obtain the last parameter cIII we first need to calculate (Δ(2iω0,τ0))-1.Let us denoteM=Δ(2iω0,τ0) and calculate the first term of this matrix; that is,(32)2iω0+e-2iω0τ0-η=2iω0+β1-iω02-η=β12-ω02-η+2iω01-β1.Hence,(33)det⁡M=β12-ω02-η+2iω01-β12-η+β2=β12-ω02+β+2iω01-β1β12-ω02-β1+2iω01-β1,M-1=1det⁡Mβ12-ω02-η+2iω01-β1η+βη+ββ12-ω02-η+2iω01-β1.Next we calculate Ψ2(h)=e2iω0hV, where(34)V=M-1d12GΦ,Φ=M-14ϵr35-3r41+r434ϵr35-3r41+r43=4ϵr35-3r4/1+r43β12-ω02+β+2iω01-β1det⁡M11=4ϵr35-3r41+r43β12-ω02-β1+2iω01-β111.Finally, we evaluate(35)d12GΨ2,Φ¯=4ϵr35-3r41+r43β12-ω02-β1+2iω01-β14ϵr35-3r41+r434ϵr35-3r41+r43=16ϵ2r65-3r421+r46β12-ω02-β1+2iω01-β111,and therefore(36)cIII=12a,a·d12GΨ2,Φ¯=16ϵ2r65-3r421+r46β12-ω02-β1+2iω01-β1a.Now, we are in a position to check the sign ofRe⁡c, where c=cI+cII+cIII. Clearly,(37)c=4ϵr21+r4435r8-22r4+5+8ϵr45-3r421+r421-β1+4ϵr45-3r42β12-ω02-β1-2iω01-β11+r42β12-ω02-β12+4ω021-β12a,and therefore(38)1+r444ϵr2Re⁡c=35r8-22r4+5+8ϵr45-3r421+r421-β1+4ϵr45-3r42β12-ω02-β11+r42β12-ω02-β12+4ω021-β12Re⁡a+8ω0ϵr45-3r421-β11+r42β12-ω02-β12+4ω021-β12Im⁡a,where Re⁡a=(1-τ0β1)/21-τ0β12+ω02τ02 and Im⁡a=-ω0τ0/21-τ0β12+ω02τ02.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 r¯<0.5, the first term is positive and of the order of units, suggesting that Re⁡c>0. 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 Re⁡c, as it is difficult to guess its sign.Now, we calculateRe⁡c for the reference parameters with ϵ=1. Let us recall that r=0.4115, β=0.02787, η=0.05418, ω0=0.9907, τ0=1.4476, and moreover β1=2η+β=0.13623. We obtain Re⁡a≈0.1486 and Im⁡a≈-0.2655. Consecutive terms in the brackets are equal to 13.11988735, 6.060143781, and -0.6953365688, while the last fraction is 1.082692850. Eventually,(39)Re⁡c≈0.604908396418.48469456·0.1486-1.082692850·0.2655≈1.487693962>0,which means that, in our reference case, the bifurcation issubcritical.In order to complement the analysis above, we can also show that in the reference case the second bifurcation appearing forτ0II is subcritical as well. We again denote τ0II by τ0, to shorten the notation. Now, the pair ω0, τ0 satisfies iω0+β+exp⁡(-iω0τ0)=0 yielding iω0+exp⁡(-iω0τ0)=-β, and from this relation we obtain the vector p=(1,-1)T. Hence, Φ(h)=exp⁡(iω0h)(1,-1)T is the eigenvector for the eigenvalue iω0. Next, we find q such that qd1Δ(iω0,τ0)p=1. As before, coordinates of q satisfies the same relation as for p, that is, q=(a,a), a==(1/2)((1+τ0β-iω0τ0)/((1+τ0β)2+ω02τ02)). Calculating the denominator of μ2 we obtain(40)d2Δiω0,τ0=-ω02+iω0β00-ω02+iω0β,next(41)q·d2Δiω0,τ0·p=2a-ω02+iω0β=iω0τ0-τ0β+1ω02τ02+τ0β+12ω02-iω0β=-ω0+iβ+τ0β2+ω02ω02τ02+τ0β+12ω0,and we easily see that the real part of this expression is negative.Next, we calculate(42)cI=12q·d13GΦ,Φ,Φ¯=12a,-a·12ϵr21+r42-12ϵr21+r42=12ϵr21+r42a,cII=-4ϵr31+r4211-β1a,-a·-4ϵr31+r424ϵr31+r42=4ϵr321+r442a1-β1,cIII=4r3ϵ2u1+r442iω0β+1+β2-ω02-β1.Eventually,(43)c=4ϵr21+r4431+r42+8ϵr41-β1+4ϵr4β2-ω02-β1-2iω0β+14ω02β+12+β2-ω02β12a,and hence,(44)sign⁡Re⁡c=sign⁡1+τ0β31+r42+8ϵr41-β1+4ϵr4β2-ω02-β14ω02β+12+β2-ω02-β12-8ϵr2τ0ω02β+14ω02β+12+β2-ω02-β12.For our reference values the expression in the brackets equals 3.146913951>0.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 I2(t)≡I and(45)I1t=I+σ,for0≤t≤Tstim,I,fort≥Tstim,where Tstim 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 1 according to the psychometric function (2) will be.In all the numerical experiments, we have chosenϵ=1 and α=3 (for this value the time scale τr 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 I=0.4 to reflect that there is a certain baseline input to the network. Other parameter values fixed across all simulations were β=100, γ=0.001, and Tstim=0.5 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 andAbsTol equal to 1e-8). 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 I.We start our numerical experiments from simulations in which stimulus of a fixed magnitudeσ=0.05 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, p1 would cross the threshold of 1-γ 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 σ=0.05. For small delays, that is, in the stability regime, we observe vanishing oscillations with no certain decision being made (p1<1). Interestingly, for τ=0.6 which is above critical destabilizing delay value we observe multiple preference switches (middle column). For τ=0.9 certain decision is being made (p1=1).Interestingly, for delays above the critical value at which first Hopf bifurcation occurs, we observe very nonintuitive behavior; see middle panels in Figures3(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 p1>p2 into p2<p1. For even larger delays we observe that a certain decision has been made for the same stimulus of σ=0.05 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 1. 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).Figure 4
(a) Decision map resulting from the psychometric function evaluation on solutions to (1) for different values of time delay τ and stimulus strengths σ. Yellow regions indicate values for which certain decision is being made in favor of population r1, that is, p1=1 after breaching critical γ distance first. Dark blue regions indicate values for which certain decision is being made in favor of population 2, that is, p2=1 after breaching critical γ distance first. Values in other places are p1(T) for T=15. (b) Number of decision switches (see a particular example in middle panel Figure 3(b)) for different delay values τ and stimulus strengths σ.
(a)
(b)In Figure5, in order to complete presentation of the numerical experiments, we show the solution for which decision made by the network is in favor of the population 2 instead of the population 1 to which the stimulus was applied.Figure 5
Exemplary solution of (1) in which stimulus is being applied to the population 1, but certain decision is being made in favor of the population 2, that is, p2=1 after breaching γ distance at some point.


## 5. Discussion

In this study, we model perceptual decision-making in elderly individuals. Although there is an extensive evidence that the cognitive impairments in aging relate to the slow processing of information [13, 16, 30, 31], the cortical mechanisms underlying cognitive impairments in aging remain elusive. In principle, the decision-making experiments in humans are mostly performed in young and healthy individuals [32, 33], which is one of the reasons behind the lack of knowledge upon the cortical mechanisms of aging.In order to bridge this gap, we propose to study the influence of synaptic delays on making perceptual choices with use of a population model of decision-making based on a winner-take-all mechanism. In this simple model, the network needs to make a binary choice between the two options. The inspiration for this model was a spiking neuronal network model with a slow reverberation mechanism by Wang [18], created to explain Shadlen and Newsome’s experiments on perceptual decision-making in rhesus monkeys [17]. Unlike Wang, we do not use simulations of stochastic spiking neuronal networks in our study though, but we take an analytic approach instead. This allows us to study the rich dynamic repertoire of networks with delay and the influence of the delay on the outcome decisions.We achieve two main results that can contribute to the understanding of the associations between the delayed GABA-signaling and the cognitive impairments during aging. Firstly, the decision-making performance is dependent on the synaptic delays in the local inhibition. For the delays below the critical value, we observe a clear association between the magnitude of the stimulus and the probability of making the correct decision (Figure4). This result is concordant with classic experiments on perceptual decision-making, as typically, a monotonic psychometric curve links the strength of the stimulus with the decision accuracy [17, 34]. However, for the synaptic delays exceeding the critical threshold, the system falls into a new dynamical regime and is no longer precise in making decisions. Dependent on both the synaptic delays and the signal magnitude, it can even achieve an accuracy of zero, by always choosing the wrong option. As known from experiments on reversal learning in rhesus monkeys [35] and humans [36], aging subjects exhibit habitual behaviors and are reluctant to relearn rules. Therefore in some cases, falling into a behavioral schema in which the wrong decision is being taken in repetitive fashion is possible in elderly subjects.The second important result from our study is that, for delay values above critical threshold, the system exhibits ambiguity in decision-making, reflected by decision switches. This result can account for the slow reaction times in perceptual decision-making in the elderly [37–41]; however the subsequent experimental validation is necessary in order to test this hypothesis.Our model predicts that impairment in the local inhibition in the cortex can result in the impaired decision-making. Although GABA concentration in prefrontal cortex and perceptual decision-making are both affected by aging, there is a lack of computational models characterizing the causal link between the two. Therefore, the model should be validated in laboratory conditions. The prediction given by the model is hard to test in the human cohorts because recordings from interneurons in the cortex are invasive. However, there are now tools in translational psychiatry that make this validation viable. For instance, the decision-making quality can be evaluated in mice in multiple experimental paradigms [19] and that results can be than correlated with the speed of synaptic transmission evaluated postmortem in an in vitro experiment [42].One remark to make is that our model is qualitative rather than quantitative, and the same mechanism can be encountered at different time scales of inhibition and for different configurations of inputs and stimuli. There are multiple GABAergic receptors in the cortex, and they have their own characteristic timescales. In example, the fast mode of inhibition is related to the GABA-A receptors [43, 44] which give synaptic delays lasting for several tens of milliseconds (excluding the afterdepolarizations lasting for several tens of milliseconds as well [45–49]). On the other hand, the slowest mode of inhibition is related to the metabotropic GABA-B receptors which have a time scale of a few hundred milliseconds [50]. Still little is known about the structure and functions of these receptors [51]. In practice, the local inhibition in the nodes of the cortical network is most probably a combination of the multiple interacting processes at different time scales. Our model is a demonstration of the principle and does not require specification of GABAergic receptors that lead a primary role in the oscillatory mechanisms considered in this work.As a summary, the model proposed in this work yields new insights into the mechanisms of aging in cortical circuits, mediated by neurodegeneration in the local inhibitory synapses.



---

*Source: 4391587-2017-09-06.xml*