$$\delta\_{n-j}^{\infty} = 0 + \varkappa \, V\_{n-j}^{\infty} - 0 = \varkappa \, V\_{n-j}^{\infty}$$

$$= a^j \varkappa^j \nu^j R / \{1 - \varkappa \,(1 - \alpha)\}^j.$$

The solid lines in **Figure 2C** show δ∞ *i* for all the states from the start (*S*1) to the goal (*S*7) when there are 7 states (*n* = 7), with varying the learning rate (α) (**Figure 2Ca**), time discount factor (γ ) (**Figure 2Cb**), decay factor () (**Figure 2Cc**), or the amount of reward (*R*) (**Figure 2Cd**) (unvaried parameters in each panel were set to the middle values: α = 0.6, γ = 0.8(1/6) , = 0.75, and *R* = 1); the dashed lines show δ∞ *i* in the model without incorporating the decay for comparison. As shown in the figures, in the cases with decay, the eventual (asymptotic) values of RPE after the convergence of learning entail gradual ramping toward the goal under a wide range of parameters. Also notably, as appeared in the " = 0.87" line in **Figure 2Cc**, depending on parameters, a peak at the start and a ramp toward the goal could coexist.

#### **MAGNITUDE-DEPENDENT RATE OF THE DECAY OF LEARNED VALUES**

We also considered cases where the rate of decay of learned values depends on the current magnitude of values so that larger values are more resistant to decay. We constructed a time-step-based model, in which decay with such magnitude-dependent rate was incorporated. Specifically, we again considered a model of the same I-maze task (**Figure 2A**) and assumed that RPE is computed at each time step *t* as follows:

$$\delta(t) = \mathcal{R}(t) + \mathcal{V}(\mathcal{S}(t)) - V(\mathcal{S}(t-1)),$$

where *S*(*t*) is the state at time step *t* and *V*(*S*(*t*)) is its value, and *R*(*t*) and δ(*t*) are obtained reward and RPE at time step *t*, respectively. γ is the time discount factor (per time step). According to this RPE, the value of state *S*(*t* − 1) was assumed to be updated as follows:

$$V\left(\mathcal{S}\left(t-1\right)\right) \to V\left(\mathcal{S}\left(t-1\right)\right) + \alpha\delta(t),$$

where α is the learning rate. We then considered the following function of value *V*:

$$\varkappa(V) = 1 - (1 - \varkappa\_1) \exp(-V/\varkappa\_2),$$

where 1 and 2 are parameters, and assumed that the value of every state decays at each time step as follows:

*V* → ( (*V*))1/*n* × *V* (for the value of states without update according to RPE), or

*V* → ( (*V*))1/*n* × (*V* + αδ) (for the value of state with update according to RPE).

**Figure 3Ba** shows the function (*V*) with a fixed value of 1 (1 = 0.6) and various values of 2 [2 = ∞ (lightest gray lines), 1.5 (second-lightest gray lines), 0.9 (dark gray lines), or 0.6 (black lines)], and **Figure 3Bb** shows the decay of learned values with each of these cases (with 7 time steps per trial assumed). For each of these cases, we simulated 100 trials of the I-maze task shown in **Figure 2A** with 7 states, with assuming γ = 0.8(1/6) and α = 0.5 and without considering reward expectation over multiple trials, and the eventual values of RPE are presented in the solid lines in **Figure 3Bc**. Notably, the time-step-based model described in the above is not exactly the same as the trial-based model described in the previous section even for the case where the rate of decay is constant: in the time-step-based model, upon the calculation of RPE: δ (*t*) = *R* (*t*) + γ*V* (*S* (*t*)) − *V* (*S* (*t* − 1)), *V* (*S* (*t*)) has suffered decay (*n* − 1) times, rather than *n* times (which correspond to a whole trial), after it has been updated last time.

#### **SIMULATION OF MAZE TASKS WITH REWARDED AND UNREWARDED GOALS**

As a simplified model of the T-maze free-choice task with rewarded and unrewarded goals used in the experiments (Howe et al., 2013) (see the Results for explanation of the task), we considered a free-choice task as illustrated in **Figure 4A**, where each state represents a relative location on the path expected to lead

**FIGURE 3 | Decay of learned values with magnitude-dependent rate leads to sigmoidal ramping of RPE resembling the observed DA ramping. (A)** was reprinted by permission from Macmillan Publishers Ltd: Nature (Howe et al., 2013), copyright (2013). **(A)** DA ramping in the ventromedial striatum observed in the experiments (Howe et al., 2013). **(B) (a)** Presumed magnitude-dependence of the rate of decay of learned values raised to the power of the number of time steps in a trial. *(Continued)*

#### **FIGURE 3 | Continued**

The horizontal black dashed line at 1 represents the case without decay, and the horizontal lightest-gray solid line at 0.6 represents the case of decay with a constant (magnitude-independent) rate. The three curved lines indicate three different degrees of magnitude-dependence of the rate of decay. **(b)** Decay of learned values under the different degrees of magnitude-dependence of the rate of decay [line colors (brightnesses) correspond to those in panel **(a)**]. **(c)** The solid lines indicate the values of RPE after 100 trials at all the states from the start (*S*1) to the goal (*S*7) in the simulated I-maze task shown in **Figure 2A** with 7 states (*n* = 7) in the model incorporating the decay with magnitude-dependent/independent rate, with varying the magnitude-dependence [line colors (brightnesses) correspond to those in **(a,b)**]. The dashed line shows the case of the model without decay.

to, or the path after passing, the rewarded or unrewarded goal or at either of the goals in each trial. We assumed that subject moves to the neighboring state in each time step, and chooses one of the two possible actions (leading to one of the two goals) at the branch point (*S*5), while learning the values of each stateaction pair (*A*1, *A*2, ··· : there is assumed to be only a single action "moving forward" in the states other than the branch point), according to one of the major reinforcement (TD) learning algorithms called Q-learning (Watkins, 1989) (for the reason why we have chosen Q-learning, see the Results), with additionally incorporating the decay of learned values with magnitudedependent rate. Specifically, at each time step *t*, RPE is computed as follows:

$$\delta(t) = R(t) + \gamma Q(A\left(t\right)) - Q\left(A\left(t - 1\right)\right) \text{(at states other than } S\_5\text{)}$$

$$\delta(t) = R(t) + \gamma \max\left\{ Q(A\_5), Q(A\_6) \right\} - Q(A\left(t - 1\right)) \text{ (at } \mathbb{S}\_5\text{)},$$

where *A*(*t*) is the state-action pair at time step *t* and *Q*(*A*(*t*)) is its value, and γ is the time discount factor (per time step). There were assumed to be *N* = 25 time steps per trial, including the inter-trial interval, and γ was set to γ = 0.81/25. According to this RPE, the value of the previous state-action pair is updated as follows:

$$Q\left(A\left(t-1\right)\right) \to Q\left(A\left(t-1\right)\right) + a\delta\left(t\right),$$

where α is the learning rate and it was set to 0.5. We then assumed that the value of every state-action pair (denoted as *Q*) decays at each time step as follows:

$$Q \to \left( \varkappa \left( Q \right) \right)^{1/N} \times Q,$$

where (*Q*) is the function introduced above, and 1 and 2 were set to 1 = 0.6 and 2 = 0.6. At the branch point (*S*5), one of the two possible actions (*A*5 and *A*6) is chosen according to the following probability:

$$\text{Prob}\left(A\_5\right) = 1/\left(1 + \exp\left(-\beta\left(Q\left(A\_5\right) - Q\left(A\_6\right)\right)\right)\right),$$

$$\text{Prob}\left(A\_6\right) = 1/\left(1 + \exp\left(-\beta\left(Q\left(A\_6\right) - Q\left(A\_5\right)\right)\right)\right)$$

$$= 1 - \text{Prob}\left(A\_5\right),$$

where Prob(*A*5) is the probability that action *A*5 is chosen, and β is a parameter determining the degree of exploration vs. exploitation upon choice (as β becomes smaller, choice becomes more and more exploratory); β was set to 1.5. In the simulations of this model, we considered reward expectation over multiple trials, specifically, we assumed that at the first time step in every trial, subject moves from the last state in the previous trial to the first state in the current trial, and RPE computation and value update are done in the same manner as in the other time steps.

In addition to the simulations of the Q-learning model, we also conducted simulations of the model with a different algorithm called SARSA (Rummery and Niranjan, 1994) (the results shown in **Figure 4F**), for which we assumed the following equation for the computation of RPE at the branch point (*S*5):

$$
\delta(t) = R(t) + \mathcal{Y}Q(A\_{\text{chosen}}) - Q(A(t-1)),
$$

where *A*chosen is the action that is actually chosen (either *A*5 or *A*6), instead of the equation for Q-learning described above. In the simulations shown in **Figure 4C**, reward *R*(*t*) was assumed to be 1 only at one of the goals (*S*8) and set to 0 otherwise, whereas in the simulations shown in **Figure 4E** and **Figure 4F**, R(*t*) was assumed to be 1 and 0.25 at the two goals (*S*8 and *S*9, respectively) and set to 0 otherwise. In addition to the modeling and simulations of the free-choice task, we also conducted simulations of a forced-choice task, which could be regarded as a simplified model of the forced-choice task examined in the experiments (Howe et al., 2013). For that, we considered sequential movements and action selection in the same state space (**Figure 4A**) but randomly determined choice (*A*5 or *A*6) at the branch point (*S*5) in each trial rather than using the choice probability function described above (while RPE of the Q-learning type, taking the max of *Q*(*A*5) and *Q*(*A*6), was still assumed), and reward *R*(*t*) at the two goals were set to 1 (large reward) and 0.25 (small reward). In each of the conditions, 1000 trials were simulated, with initial values of *Q*(*A*) set to 0 for every state-action pair *A*. We did not specifically model sessions, but we considered that the 1000 trials were divided into 25 "pseudo-sessions," each of which consists of 40 trials, so as to calculate the average and s.e.m. of the mean RPE in individual pseudo-sessions across the 25 pseudo-sessions, which are shown in the solid and dashed lines in **Figures 4Ca,Db,Ea,Fb** (in these figures, the average ± standard deviation of RPE in individual trials across trials are also shown in the error bars). **Figures 4Cb,Eb** show the RPE in 401st ∼ 440th trials. In the simulations of 1000 trials for **Figures 4C,D,E** by the Q-learning model with decay, negative RPE did not occur. By contrast, negative RPE occurred rather frequently in the SARSA model (**Figure 4F**). The ratio that the rewarded goal (*S*8) was chosen (i.e., ratio of correct trials) was 65.6, 64.5, and 64.5% in the simulations of 1000 trials for **Figures 4C,E,F**, respectively. The simulations in the present work were conducted by using MATLAB (MathWorks Inc.), and the program codes will be submitted to the ModelDB (https:// senselab.med.yale.edu/modeldb/).

#### **RESULTS**

#### **DECAY OF PLASTIC CHANGES OF SYNAPSES LEADS TO RAMPING OF RPE-REPRESENTING DA SIGNAL**

We will first show how the standard reinforcement learning algorithm called the TD learning (Sutton and Barto, 1998) works and what pattern of RPE is generated by using a virtual reward learning task, and thereafter we will consider effects of possible decay

**(A)** Simulated free-choice T-maze task with rewarded and unrewarded

*(Continued)*

#### **FIGURE 4 | Continued**

specifically modeled, and thus there does not exist a particular state that corresponds to the start of each trial. **(B)** Temporal evolution of the DA concentration in the ventromedial striatum in the experiments (Howe et al., 2013). **(a)** Average DA for rewarded (blue) or unrewarded (red) trials. **(b)** Individual trials. **(C)** Temporal evolution of the RPE in the simulations of the model incorporating the decay of learned values with magnitude-dependent rate. **(a)** The thick solid blue and red lines indicate the average, across 25 "pseudo-sessions" (see the Methods), of the mean RPE for rewarded and unrewarded trials in each pseudo-session consisting of 40 trials, respectively. The dotted lines (nearly overlapped with the solid lines) indicate these averages ± s.e.m. across pseudo-sessions. The error bars indicate the average ± standard deviation of RPE in individual trials across trials. The vertical dotted, dashed, and solid gray lines correspond to the lines in **(A)**, indicating *S*1, *S*5 (branch point), and *S*8 or *S*9 (goal) in the diagram, respectively. **(b)** Examples of the temporal evolution of RPE in individual trials in the simulations. **(Da)** DA concentration in the forced-choice task in

of plastic changes of synapses storing learned values. We considered a virtual spatial navigation task as illustrated in **Figure 2A**. In each trial, subject starts from *S*1, and moves to the neighboring state in each time step until reaching the goal (*Sn*), where reward *R* is obtained (unbranched "I-maze," rather than branched "T-maze," was considered first for simplicity). Based on the prevailing theories of neural circuit mechanisms for reinforcement learning (Montague et al., 1996; Doya, 2000), we have made the following assumptions: (1) different spatial locations, or "states," denoted as *S*1 (=start), *S*2, ··· , *Sn* (=goal, where reward *R* is obtained), are represented by different subpopulations of neurons in the subject's brain (hippocampus and/or cortical regions connecting with it), and (2) "values" of these states are stored in the changes (from the baseline) in the strength of synapses between the state-representing neurons in the cortex/hippocampus and neurons in the striatum (c.f. Pennartz et al., 2011), and thereby the value of a given state *S*, denoted as *V*(*S*), is represented by the activity of a corresponding subpopulation of striatal neurons. We have further assumed, again based on the current theories, that the following pair of computations are carried out at each state (*Si*, *i* = 1, 2,..., *n*) in the DA-CBG system: (I) DA neurons receive (indirect) impacts from the striatal neurons through basal ganglia circuits, and compute the TD RPE: δ*i* = *Ri* + γ*Vi* − *Vi* − 1, where *Ri* is reward obtained at S*i* (*Ri* = 0 unless *i* = *n*); *Vi* and *Vi* − 1 are the "values" (meaning reward expectations after leaving the states) of state *Si* and *Si* − 1, respectively; and γ (0 ≤ γ ≤ 1) is a parameter defining the degree of temporal discount of future rewards called the time discount factor, and (II) the RPE is used to update the value of the previous state (i.e., *Si* − 1) through DA-dependent plastic changes of striatal synapses: *Vi* − 1 → (*Vi* − 1 + αδ*i*), where α (0 ≤ α ≤ 1) represents the speed of learning called the learning rate, and (0 ≤ ≤ 1) is a parameter for the time-dependent decay; we first considered the case of the standard reinforcement learning model without decay (the case with = 1).

Assume that initially subject does not expect to obtain reward after completion of the maze run in individual trials and thus the "values" of all the states are 0. When reward is then introduced into the task and subject obtains reward *Rn* = *R* at the goal (*Sn*), positive RPE δ*n* = *R* + γ*Vn* − *Vn* − 1 = *R* + 0 − 0 = *R* the experiments (Howe et al., 2013). The left red vertical line indicates the branch (choice) point, while the right red line indicates another (unbranched) turning point in the M-maze used in the experiments. **(b)** RPE in the simulations of the simplified forced-choice task by the model. Configurations are the same as those in **(Ca)** except for the colors: light-green and dark-green indicate the large-reward and small-reward cases, respectively. **(E)** RPE in another set of simulations, in which it was assumed that goal-reaching (trial completion) is in itself internally rewarding, specifically, *R*(*t*) in the calculation of RPE (δ(*t*)) at the rewarded goal and the unrewarded goal was assumed to be 1 (external + internal rewards) and 0.25 (internal reward only) [rather than 1 and 0 as in the case of **(C)**], respectively. Configurations are the same as those in **(C)**. **(F) (a)** DA concentration in the dorsolateral striatum in the experiments (Howe et al., 2013). **(b)** RPE in the model incorporating the algorithm called SARSA instead of Q-learning, which was assumed in the simulations shown in **(C,Db,E)**. It was assumed that goal-reaching (trial completion) is in itself internally rewarding in the same manner as in **(E)**. Configurations are the same as those in **(Ca)**.

occurs, and it is used to update the value of *Sn* − 1 : *Vn* − 1 → 0 + αδ*n* = α*R*. Then, in the next trial, subject again obtains reward *R* at the goal (S*n*) and positive RPE occurs; this time, the RPE amounts to δ*n* = *R* + γ*Vn* − *Vn* − 1 = *R* + 0 − α*R* = (1 − α) *R*, and it is used to update the value of *Sn* − 1 : *Vn* − 1 → α*R* + αδ*n* = - 2α − α2 *R*. In this way, the value of *Sn* − 1 (*Vn* − 1) gradually increases from trial to trial, and accordingly RPE occurred at the goal (δ*n* = *R* − *Vn* − 1) gradually decreases. As long as *Vn* − 1 is smaller than *R*, positive RPE should occur and *Vn* − 1 should increase in the next trial, and eventually, *Vn* − 1 converges to *R*, and RPE (δ*n*) converges to 0 (**Figure 2Ba**) (see the Methods for mathematical details). Similarly, values of the preceding states except for the initial state (V*n* − 1, *Vn* − 2, ··· ; except for *V*1) also converge to *R* and RPE at these states (δ*n* − 1, δ*n* − 2, ··· ; except for δ1) converges to 0. Thus, from the prevailing theories of neural circuit mechanisms for reinforcement learning, it is predicted that DA neuronal response at the timing of reward and the preceding timings except for the initial timing, representing the RPE δ*n*, δ*n* − 1, δ*n* − 2, ··· , appears only transiently when reward is introduced into the task (or the amount of reward is changed), and after that transient period DA response appears only at the initial timing, as shown in the dashed lines in **Figure 2C**, which indicate eventual (asymptotic) values of RPE in the case with 7 states, with various parameters. The gradual ramping of DA signal observed in the actual reward-associated spatial navigation task (Howe et al., 2013) therefore cannot be explained by the DA=RPE hypothesis standing on the standard reinforcement (TD) learning algorithm (Niv, 2013).

Let us now assume that DA-dependent plastic changes of synaptic strengths are subject to time-dependent decay so that learned values stored in them decay with time. Let us consider a situation where *Vn* − 1 (value of *Sn* − 1) is smaller than *R* and thus positive RPE occurs at *Sn*. If there is no decay, *Vn* − 1 should be incremented exactly by the amount of this RPE multiplied by the learning rate (α) in the next trial, as seen above (**Figure 2Ba**). If there is decay, however, *Vn* − 1 should be incremented by the amount of α × RPE but simultaneously decremented by the amount of decay. By definition, RPE (δ*n* = *R* − *Vn* − 1) decreases as *Vn* − 1 increases. Therefore, if the rate (or amount) of decay is constant, *Vn* − 1 could initially increase from its initial value 0 given that the net change of *Vn* − 1 per trial (i.e., α × RPE − decay) is positive, but then the net change per trial becomes smaller and smaller as *Vn* − 1 increases, and eventually, as α × RPE becomes asymptotically equal to the amount of decay, increase of *Vn* − 1 should asymptotically terminate (**Figure 2Bb**). Even at this asymptotic limit (approximating the situation after many trials), RPE at the goal (δ*n*) remains to be positive, because it should be equal to the amount of decay divided by α. Similarly, RPE at the timings preceding reward (δ*n* − 1, δ*n* − 2, ···) also remains to be positive (see the Methods for mathematical details). The situation is thus quite different from the case without decay, in which RPE at the goal and the preceding timings except for the initial timing converges to 0 as seen above. The solid lines in **Figure 2C** show the eventual (asymptotic) values of RPE in the I-maze task (**Figure 2A**) with 7 states in the case of the model with decay, amount of which is assumed to be proportional to the current magnitude of the state value (synaptic strength) (i.e., the *rate* of decay is constant, not depending on the magnitude), with varying the learning rate (α) (**Figure 2Ca**), the time discount factor (γ ) (**Figure 2Cb**), the decay factor (κ) (**Figure 2Cc**), or the amount of reward (*R*) (**Figure 2Cd**). As shown in the figures, under a wide range of parameters, RPE entails gradual ramping toward the goal, and the ramping pattern is proportionally scaled with the amount of reward (**Figure 2Cd**).

#### **EXPLANATION OF THE OBSERVED GRADUALLY RAMPING DA SIGNAL**

As shown so far, the experimentally observed gradual ramping of DA concentration toward the goal could potentially be explained by incorporating the decay of plastic changes of synapses storing learned values into the prevailing hypothesis that the DA-CBG system implements the reinforcement learning algorithm and DA represents RPE. In the following, we will see whether and how detailed characteristics of the observed DA ramping can be explained by this account. First, the experimentally observed ramping of DA concentration in the VMS entails a nearly sigmoidal shape (**Figure 3A**) (Howe et al., 2013), whereas the pattern of RPE/DA ramping predicted from the above model (**Figure 2C**) is just convex downward, with the last part (just before the goal) being the steepest. We explored whether this discrepancy can be resolved by elaborating a model. In the model considered in the above, we assumed decay with a constant (magnitudeindependent) rate. In reality, however, the rate of decay may depend on the magnitude of learned values (synaptic strengths storing the values). Indeed, it has been shown in hippocampal slices that longer tetanus trains cause a larger degree of longterm potentiation, which tends to exhibit less decay (Gustafsson et al., 1989). Also, in the experiments examining (presumably) direct inputs from the hippocampal formation (subiculum) to the nucleus accumbens (Figure 6A of Boeijinga et al., 1993), decay of potentiation appears to be initially slow and then accelerated. We constructed an elaborate model incorporating decay with magnitude-dependent rate, which could potentially be in accord with these findings. Specifically, in the new model we assumed that larger values (stronger synapses) are more resistant to decay (see the Methods for details). We simulated the I-maze task (**Figure 2A**) with this model, and examined the eventual values of RPE after 100 trials, with systematically varying the magnitude-dependence of the rate of decay (**Figures 3Ba,b**). **Figure 3Bc** shows the results. As shown in the figure, as the magnitude-dependence of the rate of decay increases so that larger values (stronger synapses) become more and more resistant to decay, the pattern of RPE ramping changes its shape from purely convex downward to nearly sigmoidal. Therefore, the experimentally observed nearly sigmoidal DA ramping could be better explained by tuning such magnitude-dependence of the rate of decay.

Next, we examined whether the patterns of DA signal observed in the free-choice task (Howe et al., 2013), specifically, cue (tone)—reward association T-maze task, can be reproduced by our model incorporating the decay. In that task, subject started from the end of the trunk of letter "T". As the subject moved forward, a cue tone was presented. There were two different cues (1 or 8 kHz) indicating which of the two goals lead to reward in the trial. Subject was free to choose either the rewarded goal or the unrewarded goal. In the results of the experiments, subjects chose the rewarded ("correct") goal in more than a half (65%) of trials overall, indicating that they learned the cue-reward association and made advantageous choices at least to a certain extent. During the task, DA concentration in the VMS was shown to gradually ramp up, in both trials in which the rewarded goal was chosen and those in which the unrewarded goal was chosen, with higher DA concentration at late timings observed in the rewarded trials (**Figure 4Ba**). We tried to model this task by a simplified free-choice task as illustrated in **Figure 4A**, where each state represents a relative location on the path expected to lead to, or the path after passing, the rewarded or unrewarded goal or at either of the goals in each trial. The VMS, or more generally the ventral striatum receives major dopaminergic inputs from the DA neurons in the ventral tegmental area (VTA), whose activity pattern has been suggested (Roesch et al., 2007) to represent a particular form of RPE defined in one of the major reinforcement (TD) learning algorithms called Q-learning (Watkins, 1989). Therefore, we simulated sequential movements and action selection in the task shown in **Figure 4A** by using the Q-learning model incorporating the decay of learned values with magnitude-dependent rate (see the Methods for details).

Given that the model's parameters are appropriately tuned, the model's choice performance can become comparable to the experimental results (about 65% correct), and the temporal evolution of the RPE averaged across rewarded trials and also the average across unrewarded trials can entail gradual ramping during the trial (**Figure 4Ca**), reproducing a prominent feature of the experimentally observed DA signal. In the experiments (Howe et al., 2013), the authors have shown that the momentto-moment level of DA during the trial is likely to reflect the proximity to goal (location in the maze) rather than elapsed time. Although our model does not have description of absolute time and space, the value of RPE in our model is uniquely determined depending on the state, which is assumed to represent relative location in the maze, and thus given that the duration of DA's representation of RPE co-varies with the duration spent in each state, our model could potentially be consistent with the observed insensitivity to elapsed time. A major deviation of the simulated RPE/DA from the experimentally observed DA signal is that difference between rewarded trials and unrewarded trials is much larger in the simulation results, as appeared in **Figure 4Ba** and **Figure 4Ca**. We will explore how this could be addressed below. **Figure 4Cb** shows examples of the temporal evolution of RPE in individual trials in the simulations. As appeared in the figure, ramping can occur in a single trial at least for a certain fraction of trials, although more various patterns, including ramping peaked at earlier times, transient patterns, and patterns with more than one peaks, also frequently appear (see **Figure 4Bb** for comparison with the experimental results). Closely looking at the simulation results (**Figure 4Cb**), there exist oblique stripe patterns from top right to bottom left (especially clearly seen for blue colors), indicating that upward or downward deviation of RPE values, first occurred at the timing of goal and the preceding timing due to presence or absence of reward, transmits to earlier timing (to the left in the figure) in subsequent trials (to the bottom). The reason for the appearance of such a pattern is that RPE is used to update the value of state-action pair at the previous timing. This pattern is a prediction from the model and is expected to be experimentally tested, although the difference in DA signal around the timing of goal between rewarded and unrewarded trials was much smaller in the experiments, as mentioned above, and thus finding such a pattern, even if exist, would not be easy.

In the study that we modeled (Howe et al., 2013), in addition to the free-choice task, the authors also examined a forced-choice task, in which subject was pseudo-randomly forced to choose one of the goals associated with high or low reward in each trial. The authors have then found that DA ramping was strongly biased toward the goal with the larger reward (**Figure 4Da**). We considered a simplified model of the forced-choice task, represented as state transitions in the diagram shown in **Figure 4A** with the two goals associated with large and small rewards and the choice in each trial determined (pseudo-)randomly (see the Methods for details). We conducted simulations of this task by using our model with the same parameters used in the simulations of the free-choice task, and found that the model could reproduce the bias toward the goal with the larger reward (**Figure 4Db**).

#### **EXPLANATION OF FURTHER FEATURES OF THE OBSERVED DA SIGNAL**

Although our model could explain the basic features of the experimentally observed DA ramping to a certain extent, there is also a major drawback as mentioned in the above. Specifically, in our simulations of the free-choice task, gradual ramping of the mean RPE was observed in both the average across rewarded trials and the average across unrewarded trials, but there was a prominent difference between these two (**Figure 4Ca**). In particular, whereas the mean RPE for rewarded trials ramps up until subject reaches the goal, the mean RPE for unrewarded trials ramps up partway but then drops to 0 after passing the branch point. In the experiments (Howe et al., 2013), the mean RPE for rewarded trials and that for unrewarded trials did indeed differentiate later in a trial (**Figure 4Ba**), but the difference was much smaller, and the timing of differentiation was much later, than the simulation results. The discrepancy in the timing could be partially understood given that our model describes the temporal evolution of RPE, which is presumably first represented by the activity (firing rate) of DA neurons whereas the experiments measured the concentration of DA presumably released from these neurons and thus there is expected to be a time lag, as suggested from the observed difference in latencies of DA neuronal firings (Schultz et al., 1997) and DA concentration changes (Hart et al., 2014). The discrepancy in the size of the difference between rewarded and unrewarded trials, however, seems not to be explained in such a straightforward manner even partially. In the following, we would like to present a possible explanation for it.

In the simulations shown in the above, it was assumed that the unrewarded goal is literally not rewarding at all. Specifically, in our model, we assumed a positive term representing obtained reward (*R*(*t*) > 0) in the calculation of RPE (δ(*t*)) at the rewarded goal, but not at the unrewarded goal [where *R*(*t*) was set to 0]. In reality, however, it would be possible that reaching a goal (completion of a trial) is in itself internally rewarding for subjects, even if it is the unrewarded goal and no external reward is provided. In order to examine whether incorporation of the existence of such internal reward could improve the model's drawback that the difference between rewarded and unrewarded trials is too large, we conducted a new simulation in which a positive term representing obtained external or internal reward (*R*(*t*) > 0) was included in the calculation of RPE (δ(*t*)) at both the rewarded goal and the unrewarded goal, with its size four times larger in the rewarded goal [i.e., *R(t)* = 1 or 0.25 at the rewarded or unrewarded goal, respectively; this could be interpreted that external reward of 0.75 and internal reward of 0.25 are obtained at the rewarded goal whereas only internal reward of 0.25 is obtained at the unrewarded goal]. **Figure 4E** shows the results. As shown in **Figure 4Ea**, the mean RPE averaged across unrewarded trials now remains to be positive after the branch point and ramps up again toward the goal (arrowheads in the figure), and thereby the difference between rewarded and unrewarded trials has become smaller than the case without internal reward. Neural substrate of the presumed positive term (*R*(*t*)) representing internal reward is not sure, but given the suggested hierarchical reinforcement learning in the CBG circuits (Ito and Doya, 2011), such inputs might originate from a certain region in the CBG circuits that controls task execution and goal setting (in the outside of the part that is modeled in the present work).

In the study that we modeled (Howe et al., 2013), DA concentration was measured in both the VMS and the dorsolateral striatum (DLS), and there was a difference between them. In the VMS, nearly constant-rate ramping starts just after the trialonset, and rewarded and unrewarded trials differentiate only in the last period, as we have seen above (**Figure 4Ba**). In the DLS, by contrast, initial ramping looks less prominent than in the VMS, while rewarded and unrewarded trials appear to differentiate somewhat earlier than in the VMS (**Figure 4Fa**). The VMS and DLS, or more generally the ventral striatum and dorsal striatum, are suggested to receive major dopaminergic inputs from the VTA and the substantia nigra pars compacta (SNc), respectively (Ungerstedt, 1971), though things should be more complicated in reality (Björklund and Dunnett, 2007; Bromberg-Martin et al., 2010). Both VTA and SNc DA neurons have been shown to represent RPE, but they may represent different forms of RPE used for different reinforcement (TD) learning algorithms. Specifically, it has been empirically suggested, albeit in different species, that VTA and SNc DA neurons represent RPE for Q-learning (Roesch et al., 2007) and SARSA (Morris et al., 2006), respectively; these two algorithms differ in whether the maximum value of all the choice options (Q-learning) or the value of actually chosen option (SARSA) is used for the calculation of RPE [see (Niv et al., 2006) and the Methods]. Conforming to this suggested distinction, so far we have assumed Q-learning in the model and compared the simulation results with the DA concentration in the VMS that receives major inputs from the VTA. The emerging question, then, is whether simulation results become more comparable to the DA concentration in the DLS if we instead assume SARSA in the model. We explored this possibility by conducting a new simulation, and found that it would indeed be the case. **Figure 4Fb** shows the simulation results of the model with SARSA, which also incorporated the internal reward upon reaching the unrewarded goal introduced above. Compared with the results with Q-leaning (**Figure 4Ea**), initial ramping looks less prominent, and rewarded and unrewarded trials differentiate earlier. These two differences could be said to be in line with the experimentally observed differences between the VMS and DLS DA concentrations as described above, although again the difference between rewarded and unrewarded trials is larger, and the timing of differentiation is earlier, in the model than in the experiment.

Intriguingly, in the study that has shown the representation of RPE for Q-learning in VTA DA neurons (Roesch et al., 2007), DA neurons increased their activity in a staggered manner from the beginning of a trial (before cue presentation) toward reward, with the activity in the middle of the increase shown to entail the characteristics of RPE. It is tempting to guess that such a staggered increase of VTA DA neuronal firing actually has the same mechanistic origin as the gradual increase of VMS DA concentration in the study that we modeled (Howe et al., 2013). Consistent with this possibility, in a recent study that has simulated the experiments in which VTA DA neurons were recorded (Roesch et al., 2007) by using a neural circuit model of the DA-CBG system (Morita et al., 2013), the authors have incorporated decay of learned values, in a similar manner to the present work, in order to reproduce the observed temporal pattern of DA neuronal firing, in particular, the within-trial increase toward reward (although it was not the main focus of that study and also the present work does not rely on the specific circuit structure/mechanism for RPE computation proposed in that study).

## **DISCUSSION**

While the hypothesis that DA represents RPE and DA-dependent synaptic plasticity implements update of reward expectations according to RPE has become widely appreciated, recent work has revealed the existence of gradually ramping DA signal that appears not to represent RPE. We explored whether such DA ramping can be explained by extending the "DA=RPE" hypothesis by taking into account possible time-dependent decay of DA-dependent plastic changes of synapses storing learned values. Through simulations of reward learning tasks by the RPE-based reinforcement learning model, we have shown that incorporation of the decay of learned values can indeed cause gradual ramping of RPE and could thus potentially explain the observed DA ramping. In the following, we discuss limitations of the present work, comparisons and relations with other studies, and functional implications.

#### **LIMITATIONS OF THE PRESENT WORK**

In the study that has found the ramping DA signal (Howe et al., 2013), it was shown that the peak of the ramping signals was nearly as large as the peak of transient responses to unpredicted reward. By contrast, in our simulations shown in **Figure 4E**, average RPE for all the trials at state *S*5 is about 0.158, which is smaller than RPE for unpredicted reward of the same size in our model (it is 1.0). This appears to deviate from the results of the experiments. However, there are at least three potential reasons that could explain the discrepancy between the experiments and our modeling results, as we describe below.

First, in the experiments, whereas there was only a small difference between the peak of DA response to free reward and the peak of DA ramping during the maze task when averaged across sessions, the slope of the regression line between these two values (DA ramping / DA to free reward) in individual sessions (Extended Data Figure 5a of Howe et al., 2013) is much smaller than 1 (it is about 0.26). Indeed, that figure shows that there were rather many sessions in which the peak of DA response to free reward was fairly large (>15 nM) whereas the peak of DA ramping during the maze task was not large (<15 nM), while much less sessions exhibited the opposite pattern. How the large variability in DA responses in the experiments reflects heterogeneity of DA cells and/or other factors is not sure, but it might be possible to regard our model as a model of cells or conditions in which response to free reward was fairly large whereas ramping during the maze task was not large. Second, it is described in Howe et al. (2013) (legend of Extended Data Figure 5a) that DA response to free reward was compared with DA ramping measured from the same probes during *preceding* behavioral training in the maze. Given that the same type of reward (chocolate milk) was used in the task and as free reward, and that the measurements of DA response to deliveries of free reward were made after the measurements of DA ramping during 40 maze-task trials in individual sessions, we would think that there possibly existed effects of satiety. Third, the degree of unpredictability of the "unexpected reward" in the experiments could matter. Specifically, it seems possible that there were some sensory stimuli that immediately preceded reward delivery and informed the subjects of it such as sounds (generated in the device for reward supply) or smells. In such a case, conventional RPE models without decay predict that, after some experience of free reward, RPE of nearly the same size as that of RPE generated upon receiving ultimately unpredictable reward is generated at the timing of the sensory stimuli (unless time discount is extremely severe: size becomes smaller only due to time discount), and no RPE is generated at the timing of actual reward delivery. In contrast, and crucially, our model with decay predicts that, after some experience of free reward, RPE generated at the timing of the sensory stimuli is significantly smaller than RPE generated upon receiving ultimately unpredictable reward, and positive RPE also occurs upon receiving reward but it is also smaller than the ultimately unpredictable case (if the timing of the sensory stimuli is one time-step before the timing of reward in the model with the parameters used for **Figure 4**, RPE values at those two timings after 15 experiences are about 0.87 and 0.16, respectively; these two RPEs are about 0.99 and 0.00 in the case without decay). The mechanism of this can be schematically understood from **Figure 2Bb** by viewing *Vn* − 1 (bar height) and δ*n* (space above the bar) as RPEs at the timings of the preceding sensory stimuli and the actual reward delivery, respectively (as for the former, except for time discount); they are both smaller than the reward amount ("*R*"), which is the size of RPE generated upon receiving this reward ultimately unpredictably. With these considerations, we would think that the discrepancy between the experiments and the model in the relative sizes of the peak DA response to free reward and the peak DA ramping in the maze task could potentially be explained.

Other than the point described above, there are at least six fundamental limitations of our model. First, our model's behavior is sensitive to the magnitude of rewards. As shown in the Results, in our original model assuming decay with a constant rate, overall temporal evolution of RPE is proportionally scaled according to the amount of reward (**Figure 2Cd**). However, such a scalability no longer holds for the elaborated model incorporating the magnitude-dependent rate of decay, because the assumed magnitude-dependence (**Figure 3Ba**) is sensitive to absolute reward amount. Consequently, the patterns of RPE shown in **Figures 3** and **4** will change if absolute magnitude of rewards is changed. In reality, it is possible that magnitude-dependence of the rate of decay of learned values (synaptic strength) itself can be changed, in a longer time scale, depending on the average magnitude of rewards obtained in the current context. Second, whereas the free-choice task used in the experiments (Howe et al., 2013) involved cue-reward association, our simplified model does not describe it. Because of this, the state in our model is assumed to represent relative location on the path expected to lead to, or the path after passing, the rewarded or unrewarded goal or at either of the goals in each trial (as described before), but not absolute location since the absolute location of rewarded/unrewarded goal in the experiments was determined by the cue, which changed from trial to trial. Third, our model only has abstract representation of relative time and space, and how they are linked with absolute time and space is not defined. Fourth, validity of our key assumption that plastic changes of synapses are subject to time-dependent decay remains to be proven. There have been several empirical suggestions for the (rise and) decay of synaptic potentiation (Gustafsson et al., 1989; Boeijinga et al., 1993; Xiao et al., 1996) and spine enlargement (Matsuzaki et al., 2004) in the time scale of minutes, which could potentially fit the time scale of the maze task simulated in the present study, but we are currently unaware of any reported evidence for (or against) the occurrence of decay of DA-dependent plastic changes of synapses in animals engaged in tasks like the one simulated in the present study. Also, we assumed simple equations for the decay, but they would need to be revised in future works. For example, any plastic changes will eventually decay back to 0 according to the models in the present work, but in reality at least some portion of the changes is likely to persist for a long term as shown in the experiments referred to in the above. Fifth, regarding the origin of the ramping DA signal and its potential relationships with the DA=RPE hypothesis, there are potentially many possibilities, and the mechanism based on the decay of learned values proposed in the present study is no more than one of them (see the next section for two of other possibilities). Sixth, potential modulation of DA release apart from DA neuronal firing is not considered in the present study. We have assumed that the observed ramping DA signal in the striatum (Howe et al., 2013) faithfully reflects DA neuronal firing, which has been suggested to represent RPE. However, as pointed out previously (Howe et al., 2013; Niv, 2013), whether it indeed holds or not is yet to be determined, because DA neuronal activity was not measured in that study and DA concentration can be affected by presynaptic modulations of DA release, including the one through activation of nicotinic receptors on DA neuronal axons by cholinergic interneurons (Threlfell et al., 2012), and/or saturation of DA reuptake. Addressing these limitations would be interesting topics for future research.

#### **COMPARISONS AND RELATIONS WITH OTHER STUDIES**

Regarding potential relationships between the ramping DA signal in the spatial navigation task and the DA=RPE hypothesis, a recent theoretical study (Gershman, 2014) has shown that DA ramping can be explained in terms of RPE given nonlinear representation of space. This is an interesting possibility, and it is entirely different from our present proposal. The author has argued that his model is consistent with important features of the observed DA ramping, including the dependence on the amount of reward and the insensitivity to time until the goal is reached. Both of these features could also potentially be consistent with our model, although there are issues regarding the sensitivity of model's behavior to reward magnitude and the lack of representation of absolute time and space, as we have so far described. It remains to be seen whether the limitations of our model, including the large difference between rewarded and unrewarded trials, are not the case with his model. Notably, these two models are not mutually exclusive, and it is possible that the observed DA ramping is a product of multiple factors. Also, the possible correspondence between the differential DA signal in the ventral vs. dorsal striatum and Q-learning vs. SARSA mentioned in the Results could also hold with Gershman's model.

It has also been shown (Niv et al., 2005) that the conventional reinforcement learning model (without decay) can potentially explain ramping of averaged DA neuronal activity observed in a task with probabilistic rewards (Fiorillo et al., 2003), if it is assumed that positive and negative RPEs are asymmetrically represented by increase and decrease of DA neuronal activity from the baseline, with the dynamic range of the decrease narrower than that of the increase due to the lowness of the baseline firing rate. This mechanism did not contribute to the ramping of RPE in our simulations, because such asymmetrical representation was not incorporated into our model; actually, negative RPE did not occur in the 1000-trials simulations of our Q-learning model for **Figures 4C,Db,E**, while negative RPE occurred rather frequently in the SARSA model (**Figure 4Fb**). Notably, according to the mechanism based on the asymmetrical RPE representation by DA (Niv et al., 2005), ramping would not appear in the I-maze task where reward is obtained in every trial without uncertainty (**Figure 2A**) because negative RPE would not occur in such a situation, different from the cases of the decay-based mechanism proposed in the present work and the mechanism proposed by Gershman (Gershman, 2014) mentioned above. Experimental examination of the I-maze would thus be potentially useful to distinguish mechanisms that actually operate. In the meantime, the mechanism based on the asymmetrical RPE representation by DA is not mutually exclusive with the other two, and two or three mechanisms might simultaneously operate in reality.

#### **FUNCTIONAL IMPLICATIONS**

Given that the observed DA ramping is indicative of decay of learned values as we have proposed, what is the functional advantage of such decay? Decay would naturally lead to forgetting, which is rather disadvantageous in many cases. However, forgetting can instead be useful in certain situations, in particular, where environments are dynamically changing and subjects should continually overwrite old memories with new ones. Indeed, it has recently been proposed that decay of plastic changes of synapses might be used for active forgetting (Hardt et al., 2013, 2014). Inspired by this, here we propose a possible functional advantage of synaptic decay specifically for the DA-CBG system involved in value learning. In value learning, active forgetting is required when associations between rewards and preceding sensory stimuli are changed, such as the case of reversal learning in which cue-reward association is reversed unpredictably. In theory, flexible reversal of leaned association should be possible based solely on RPE without any decay: old association can be erased by negative RPE first, and new association can then be learned by positive RPE. However, in reality there would be a problem due to a biological constraint. Specifically, it has been indicated that the dynamic range of DA neuronal activity toward the negative direction from the baseline firing rate is much narrower than the positive side, presumably for the sake of minimizing energy cost (c.f., Laughlin, 2001; Bolam and Pissadaki, 2012; Pissadaki and Bolam, 2013), and thereby DA neurons can well represent positive RPE, but perhaps not negative RPE (Bayer and Glimcher, 2005) (see also Potjans et al., 2011). This indication has been challenged by subsequent studies: it has been shown (Bayer et al., 2007) that negative RPE was correlated with the duration of pause of DA neuronal firing, and a recent study using FSCV (Hart et al., 2014) has shown that DA concentration in the striatum in fact symmetrically encoded positive and negative RPE in the range tested in that study. Nevertheless, it could still be possible that representation of negative RPE by DA is limited in case the baseline DA concentration is low. In such a case, synaptic decay could be an alternative or additional mechanism for erasing old, already irrelevant cue-reward associations so as to enable flexible reversal/reconstruction of associations, with possibly the rate of decay itself changing appropriately (i.e., speeding up just after the reversal/changes in the environments) through certain mechanisms (e.g., monitoring of the rate of reward acquisition). We thus propose that decay of learned values stored in the DA-dependent plastic changes of CBG (corticostriatal) synapses would be a feature of the DA-CBG circuits, which endows the reinforcement learning system with flexibility, in a way that is also compatible with the minimization of energy cost.

With such consideration, it is suggestive that DA ramping was observed in the study using the spatial navigation task (Howe et al., 2013) but not in many other studies (though there could be symptoms as we discussed above). Presumably, it reflects that the spatial navigation task is ecologically more relevant, for rats, than many other laboratory tasks. In the wild, rats navigate to forage in dynamically changing environments, where flexibility of learning would be pivotal. Moreover, the overall rate of rewards in wild foraging would be lower than in many laboratory tasks, and given the suggestion that the rate of rewards is represented by the background concentration of DA (termed tonic DA) (Niv et al., 2007), tonic DA in foraging rats is expected to be low and thus representation of negative RPE by DA could be limited as discussed above. The rate of decay of learned values would therefore be adaptively set to be high so as to turn on the alternative mechanism for flexible learning, and it would manifest as the prominent ramping of DA/RPE in the task mimicking foraging navigation (even if the rate of rewards is not that low in the task, different from real foraging). If this conjecture is true, changing the volatility of the task, mimicking changes in the volatility of the environment, may induce adaptive changes in the rate of decay of learned values (synaptic strengths), which could cause changes in the property of DA ramping (c.f., **Figure 2Cc**): a testable prediction of our model.

Apart from the decay, DA ramping can also have more direct functional meanings. Along with its roles in plasticity induction, DA also has significant modulatory effects on the responsiveness of recipient neurons. In particular, DA is known to modulate the activity of the two types of striatal projection neurons to the opposite directions (Gerfen and Surmeier, 2011). Then, given that DA neurons compute RPE based on value-representing BG inputs, on which the activity of striatal neurons have direct and/or indirect impacts, ramping DA, presumably representing a gradual increase of RPE according to our model, would modulate the activity of striatal neurons and thereby eventually affect the computation of RPE itself. Such a closed-loop effects (c.f., **Figure 1A**) can potentially cause rich nonlinear phenomena through recurrent iterations. Exactly what happens depends on the precise mechanism of RPE computation, while the present work does not assume specific mechanism for it so that the results presented so far can generally hold. Just as an example, however, when the model of the present study is developed into a model of the DA-CBG circuit based on a recently proposed mechanism for RPE computation (Morita et al., 2012, 2013; Morita, 2014), consideration of the effects of DA on the responsiveness of striatal projection neurons can lead to an increase in the ratio of correct trials, indicating occurrence of positive feedback (unpublished observation). This could potentially represent self-enhancement of internal value or motivation (c.f., Niv et al., 2007). Such an exciting possibility is also expected to be explored in future work.

### **AUTHOR CONTRIBUTIONS**

Kenji Morita conceived and designed the research. Kenji Morita and Ayaka Kato performed the modeling, calculations, and simulations. Kenji Morita drafted the manuscript. Ayaka Kato commented on the manuscript, and contributed to its revision and elaboration.

## **ACKNOWLEDGMENTS**

This work was supported by Grant-in-Aid for Scientific Research on Innovative Areas "Mesoscopic Neurocircuitry" (No.25115709) of The Ministry of Education, Science, Sports and Culture of Japan and Strategic Japanese - German Cooperative Programme on "Computational Neuroscience" (project title: neural circuit mechanisms of reinforcement learning) of Japan Science and Technology Agency to Kenji Morita.

## **REFERENCES**

Xiao, M. Y., Niu, Y. P., and Wigström, H. (1996). Activity-dependent decay of early LTP revealed by dual EPSP recording in hippocampal slices from young rats. *Eur. J. Neurosci.* 8, 1916–1923. doi: 10.1111/j.1460-9568.1996.tb01335.x

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 17 January 2014; accepted: 24 March 2014; published online: 09 April 2014. Citation: Morita K and Kato A (2014) Striatal dopamine ramping may indicate flexible reinforcement learning with forgetting in the cortico-basal ganglia circuits. Front. Neural Circuits 8:36. doi: 10.3389/fncir.2014.00036*

*This article was submitted to the journal Frontiers in Neural Circuits.*

*Copyright © 2014 Morita and Kato. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Neuromodulatory adaptive combination of correlation-based learning in cerebellum and reward-based learning in basal ganglia for goal-directed behavior control

## *Sakyasingha Dasgupta1,2\*, Florentin Wörgötter 1,2 and Poramate Manoonpong2,3*

*1 Institute for Physics - Biophysics, George-August-University, Göttingen, Germany*

*2 Bernstein Center for Computational Neuroscience, George-August-University, Göttingen, Germany*

*3 Center for Biorobotics, Maersk Mc-Kinney Møller Institute, University of Southern Denmark, Odense, Denmark*

#### *Edited by:*

*M. Victoria Puig, Massachusetts Institute of Technology, USA*

#### *Reviewed by:*

*Kenji Morita, The University of Tokyo, Japan Bernd Porr, University of Glasgow, UK*

#### *\*Correspondence:*

*Sakyasingha Dasgupta, Bernstein Center for Computational Neuroscience, George-August-University, Friedrich-Hund Platz 1, 37077 Göttingen, Germany e-mail<EMAIL_ADDRESS>
Goal-directed decision making in biological systems is broadly based on associations between conditional and unconditional stimuli. This can be further classified as classical conditioning (correlation-based learning) and operant conditioning (reward-based learning). A number of computational and experimental studies have well established the role of the basal ganglia in reward-based learning, where as the cerebellum plays an important role in developing specific conditioned responses. Although viewed as distinct learning systems, recent animal experiments point toward their complementary role in behavioral learning, and also show the existence of substantial two-way communication between these two brain structures. Based on this notion of co-operative learning, in this paper we hypothesize that the basal ganglia and cerebellar learning systems work in parallel and interact with each other. We envision that such an interaction is influenced by reward modulated heterosynaptic plasticity (RMHP) rule at the thalamus, guiding the overall goal directed behavior. Using a recurrent neural network actor-critic model of the basal ganglia and a feed-forward correlation-based learning model of the cerebellum, we demonstrate that the RMHP rule can effectively balance the outcomes of the two learning systems. This is tested using simulated environments of increasing complexity with a four-wheeled robot in a foraging task in both static and dynamic configurations. Although modeled with a simplified level of biological abstraction, we clearly demonstrate that such a RMHP induced combinatorial learning mechanism, leads to stabler and faster learning of goal-directed behaviors, in comparison to the individual systems. Thus, in this paper we provide a computational model for adaptive combination of the basal ganglia and cerebellum learning systems by way of neuromodulated plasticity for goal-directed decision making in biological and bio-mimetic organisms.

**Keywords: decision making, recurrent neural networks, basal ganglia, cerebellum, operant conditioning, classical conditioning, neuromodulation, correlation learning**

## **1. INTRODUCTION**

Associative learning by way of conditioning, forms the main behavioral paradigm that drives goal-directed decision making in biological organisms. Typically, this can be further classified into two classes, namely, classical conditioning (or correlation-based learning) (Pavlov, 1927) and operant conditioning (or reinforcement learning) (Skinner, 1938). In general, classical conditioning is driven by associations between an early occurring conditional stimulus (CS) and a late occurring unconditional stimulus (US), which lead to conditioned responses (CR) or unconditioned responses (UR) in the organism (Clark and Squire, 1998; Freeman and Steinmetz, 2011). The CS here acts as a predictor signal such that, after repeated pairing of the two stimuli, the behavior of the organism is driven by the CR (adaptive reflex action) at the occurrence of the predictive CS, much before the US arrives. The overall behavior is guided on the sole basis of stimulus-response (S-R) associations or correlations, without any explicit feedback in the form of rewards or punishments from the environment. In contrast to such classically conditioned reflexive behavior acquisition, operant conditioning provides an organism with adaptive control over the environment with the help of explicit positive or negative reinforcements (evaluative feedback) given for corresponding actions. Over sufficient time, this enables the organism to respond with good behaviors, while avoiding bad or negative behaviors. As such within the computational learning framework, this is usually termed reinforcement learning (RL) (Sutton and Barto, 1998).

At a behavioral level, although the two conditioning paradigms of associative learning appear to be distinct from each other, they seem to occur in combination as suggested from several animal behavioral studies (Rescorla and Solomon, 1967; Dayan and Balleine, 2002; Barnard, 2004). Behavioral studies with rabbits (Lovibond, 1983) demonstrate that the strength of operant responses can be influenced by simultaneous presentation of classically conditioned stimuli. This was further elaborated upon in the behavior of fruit flies (Drosophila), where both classical and operant conditioning predictors influence the behavior at the same time and in turn improve the learned responses (Brembs and Heisenberg, 2000). On a neuronal level, this relates to the interaction between the reward modulated action selection at the basal ganglia and the correlation based delay conditioning at the cerebellum. Although the classical notion has been to regard the basal ganglia and the cerebellum to be primarily responsible for motor control, increasing evidence points toward their role in non-motor specific cognitive tasks like goal-directed decision making (Middleton and Strick, 1994; Doya, 1999). Interestingly, recent experimental studies (Neychev et al., 2008; Bostan et al., 2010) show that the the basal ganglia and cerebellum not only form multi-synaptic loops with the cerebral cortex, but, two-way communication between the structures exist via the thalamus **Figure 1A**) along with substantial disynaptic projections to the cerebellar cortex from the subthalamic nucleus (STN) of the basal ganglia and from the dentate nucleus (cerebellar output stage) to the striatum (basal ganglia input stage) (Hoshi et al., 2005). This suggests that the two structures are not separate performing distinct functional operations (Doya, 2000a), but are linked together forming an integrated functional network. Such integrated behavior is further illustrated in the timing and error prediction studies of Dreher and Grafman (2002) showing that the activation of the cerebellum and basal ganglia are not specific to switching attention, as previously believed, because both these regions were activated during switching between tasks as well as during the simultaneous maintenance of two tasks.

Based on these compelling evidences we formulate the *neural combined learning hypothesis*, which proposes that goal-directed decision making occurs with a parallel adaptive combination (balancing) of the two learning systems (**Figure 1B**) to guide the final action selection. As evident from experimental studies (Haber and Calzavara, 2009), the thalamus potentially plays a critical role in integrating the neural signals from the two subnetworks while having the ability to modulate behavior through dopaminergic projections from the ventral tagmental area (VTA)

conditioned reflexive behaviors. Adapted and modified from Doya (2000a). **(B)** Combinatorial learning framework with parallel combination of ICO

(García-Cabezas et al., 2007; Varela, 2014). The motor thalamic (Mthal) relay nuclei, specifically the VA-VL (ventral anterior and ventral lateral) regions receive projections from the basal ganglia (inputs from the globas pallidus) as well as the cerebellum (inputs from the dentate nucleus) (Jones et al., 1985; Percheron et al., 1996). This can be further segregated with the ventral anterior and the anterior region of the ventrolateral nucleus (VLa) receiving major inputs from the globus pallidus internus (GPi), while the posterior region of the ventrolateral nucleus (VLp) receives primary inputs from the cerebellum (Bosch-Bouju et al., 2013). Recent studies using molecular markers were able to distinguish the VA and VL nuclei in rats (Kuramoto et al., 2009), which had hitherto been difficult and were considered as a single overlapping area as the VA-VL complex. Interestingly, despite apparent anatomical segregation of information in the basal ganglia and cerebellar territories, similar ranges of firing rate and movement related activity are observed in the Mthal neurons across all regions (Anderson and Turner, 1991). Furthermore, some experimental studies based on triple labeling techniques found zones of overlapping projections, as well as interdigitating foci of pallidal and cerebellar labels, particularly in border regions of the VLa (Sakai et al., 2000). In light of these evidences, it is plausible that the basal ganglia and cerebellar circuitries not only form an integrated functional network, but their individual outputs are combined together by a subset of the VLa neurons which in turn project to the supplementary and pre-supplementary motor cortical areas (Akkal et al., 2007) responsible for goaldirected movements. We envision that such a combined learning mechanism may be driven by reward modulated heterosynaptic plasticity (neuromodulation by way of dopaminergic projections) at the thalamus.

In this study, input correlation learning (ICO)in the form of a differential Hebbian learner (Porr and Wörgötter, 2006), was implemented as an example of delay conditioning in the cerebellum, while a reservoir network (Jaeger and Haas, 2004) based continuous actor-critic learner (Doya, 2000b) was implemented as an example of reward based conditioning in the basal ganglia. Taking advantage of the individual learning mechanisms, the combined framework can learn the appropriate goal-directed control policy

*Ocom* controls the agent behavior (policy) while sensory feedback from the

agent is sent back to both the learning mechanisms in parallel.

for an agent1 in a fast and robust manner outperforming the singular implementation of the individual components.

Although there have been a number of studies which have applied the two different conditioning concepts for studying selforganizing behavior in artificial agents and robots, they have mostly been applied separately to generate specific goal-directed behaviors (Morimoto and Doya, 2001; Verschure and Mintz, 2001; Hofstoetter et al., 2002; Prescott et al., 2006; Manoonpong et al., 2007; Soltoggio et al., 2013). In our previous work (Manoonpong et al., 2013) we motivated a combined approach of the two learning concepts on a purely algorithmic level without any adaptive combination between the two. To the best of our knowledge, in this paper we present for the first time a biologically plausible approach to model an adaptive combination of the cerebellar and basal ganglia learning systems, where they indirectly interact through sensory feedback. In this manner they work as a single functional unit to guide the behavior of artificial agents. We test our neural combined learning hypothesis within the framework of goal-directed decision making using a simulated wheeled robot situated in environments of increasing complexity designed as part of static and dynamic foraging tasks (Sul et al., 2011). Our results clearly show that the proposed mechanism enables the artificial agent to successfully learn the task in the different environments with changing levels of interaction between the two learning systems. Although we take a simplified approach of simulated robot based goal-directed learning, we believe our model covers a reasonable level of biological abstraction that can help us understand better, the closed-loop interactions between these two neural subsystems as evident from experimental studies and also provide a computational model of such combined learning behavior which has hitherto been missing.

We now give a brief introduction to the neural substrates of the cerebellum and the basal ganglia with regards to classical and operant conditioning. Using a broad high-level view of the anatomical connections of these two brain structures, we motivate how goal-directed behavior is influenced by the respective

1Agent here refers to any artificial or biological organism situated in a given environment.

structures and their associated neuronal connections. The individual computational models with implementation details of the two interacting learning systems are then presented in the Materials and Methods Section followed by results and discussion.

#### **1.1. CLASSICAL CONDITIONING IN THE CEREBELLUM**

The role of the Cerebellum and its associated circuitry in the acquisition and retention of anticipatory responses (sensory predictions) with Pavlovian delay conditioning has been well established (Christian and Thompson, 2003; Thompson and Steinmetz, 2009). Although most of the classical conditioning studies are primarily based on eye-blink conditioning (Yeo and Hesslow, 1998), recent experimental studies have established the essential role of the cerebellum in learning and memory of goal-directed behavioral responses (Burguiere et al., 2010). In **Figure 2A** a highly simplified control structure of the major cerebellar pathways and their relative function is indicated. The Inferior Olive relays the US signal to the cerebellar cortex through the climbing fibers and then induces plasticity at the synaptic junctions of the mossy fibers carrying the CS information (Herreros and Verschure, 2013). Repeated CS-US pairings gradually lead (through synaptic consolidation) to the acquisition of the CR with a drop in the firing activity of the Purkinje cells (output from the cerebellar cortex). The cerebral cortex projects to the lateral cerebellum via pontine nuclei relays (Allen and Tsukahara, 1974; Lisberger and Thach, 2013; Proville et al., 2014) which in turn have projections back to the cerebral cortex through relays in the thalamus (ventro-lateral nucleus), thus projecting the conditioned responses from the cerebellum to the motor cortical areas (Stepniewska et al., 1994; Sakai et al., 2000). In essence, the cerebellar action modulates or controls the motor activity of the animal which produces changes in its goal oriented behavior. The goal oriented behaviors can typically involve both attraction toward or avoidance of specific actions (generally referred to as adaptive reflexes) involving both sensory predictions and motor control, toward which the cerebellum makes a major contribution. It is also important to note that although numerous experimental and computational studies demonstrate

the function of the Cerebellum in classical conditioning or correlation learning (Kim and Thompson, 1997; Woodruff-Pak and Disterhoft, 2008), a possible role of the Cerebellum toward supervised learning (SL) has also been widely suggested (Doya, 1999; Kawato et al., 2011). Typically within the paradigm of SL a training or instructive signal acts as a reference toward which the output of a network (movements) is compared, such that the error generated acts as the driver signal to induce plasticity within the network in order to find the correct mapping between the sensory input stimuli and the desired outputs (Knudsen, 1994). Using the classical conditioning paradigm, it has been suggested that the instructive signal that supervises the learning is the input activity associated with the US. As such, the SL model of the cerebellum considers that the climbing fibers from the inferior olive provide the error signal (instructive activity) for the Purkinje cells. Coincident inputs from the inferior olive and the granule cells lead to plasticity at the granule-to-Purkinje synapses (Doya, 2000a). Although there have been experimental studies to validate the SL description of the cerebellum (Kitazawa et al., 1998), it has been largely directed toward considering the cerebellum as an internal model of the body and the environment (Kawato, 1999). Furthermore, Krupa et al. (1993) observed that even when the red nucleus (relay between motor cortex and cerebellum) was inactivated learning proceeded with no CR being expressed. Thus, this demonstrates that no error signal based on the behavior was needed for learning to occur. Instead, the powerful climbing fiber activity evoked by the US, acting as a template, could cause the connection strengths of sensory inputs that are consistently correlated with it to increase. Subsequently , after sufficient repetition, the activity of these sensory inputs alone would drive the UR pathway. As such, in this work we directly consider correlation learning as the basis of classical conditioning in the cerebellum without taking into consideration SL mechanisms and do not explicitly consider the US relay from the inferior olive as an error signal.

#### **1.2. REWARD LEARNING IN THE BASAL GANGLIA**

In contrast to the role of the cerebellum in classical conditioning, the basal ganglia and its associated circuitry possess the necessary anatomical features (neural substrates) required for a rewardbased learning mechanism (Schultz and Dickinson, 2000). In **Figure 2B** we depict the main anatomical connections of the cortical basal ganglia circuitry. It is comprised of the striatum (consisting of most of the caudate and the putamen, and of the nucleus accumbens), the internal (medial) and external (lateral) segments of the globus pallidus (GPi and GPe respectively), the subthalamic nucleus (STN), the ventral tegmental area (VTA) and the substantia nigra pars compacta (SNc) and pars reticulata (SNr). The input stage of the basal ganglia is the striatum connected via direct cortical projections. Previous studies have not only recognized the striatum as a critical structure in the learning of stimulusresponse behaviors, but also established it as the major location which projects to as well as receives efferent connections from (via direct and indirect multi-synaptic pathways) the dopaminergic system (Joel and Weiner, 2000; Kreitzer and Malenka, 2008). The processing of rewarding stimuli is primarily modulated by the dopamine neurons (DA system in **Figure 2B**) of the VTA and SNc with numerous experimental studies (Schultz and Dickinson, 2000) demonstrating, that changes in dopamine neurons encode the prediction error in appetitive learning scenarios, and associative learning in general (Puig and Mille, 2012). **Figure 2B**—right shows the idealized reciprocal architecture of the striatal and dopaminergic circuitry. Here sensory stimuli arrive as input from the cortex to the striatal network. Excitatory as well as inhibitory synapses project from the striatum to the DA system which in turn uses the changes in the activity of DA neurons to modulate the activity in the striatum. Such DA activity also acts as the neuromodulatory signal to the thalamus which receives indirect connections from the striatum, via the GPi, SNr and VTA (Varela, 2014). Computational modeling of such dopamine modulated reward learning behavior is particularly well reflected by the Temporal Difference (TD) algorithm (Sutton, 1988; Suri and Schultz, 2001), as well as in the action selection based computational models of the basal ganglia (Gurney et al., 2001; Humphries et al., 2006). In the context of basal ganglia modeling, Actor-Critic models (explained further in the next section) of TD learning (Houk et al., 1995; Joel et al., 2002) have been extensively used. They create a functional separation between two sub-networks of the critic (modeling striatal and dopaminergic activity) and the actor (modeling striatal to motor thalamus projections). The TD learning rule uses the prediction error (TD error) between two subsequent predictions of the net weighted sum of future rewards based on current input and actions, to modulate critic weights via long-term synaptic plasticity. The same prediction error signal (dopaminergic projections) is also used to modulate the synaptic weights at the actor; output from which controls the the actions taken by the agent. Typically, here the mechanism of action selection, can be regarded as the neuromodulation process occurring at the striatum, which then reaches the motor thalamic regions via projections from the output stages of the basal ganglia, namely GPi/GPe and SNr (Gurney et al., 2001; Houk et al., 2007) (**Figure 2B**).

## **2. MATERIALS AND METHODS**

### **2.1. COMBINATORIAL LEARNING WITH REWARD MODULATED HETEROSYNAPTIC PLASTICITY**

According to the neural combined learning hypothesis for successful goal-directed decision making, the underlying neural machinery of animals combines basal ganglia and cerebellar learning systems output, induced with a reward modulated balancing (neuromodulation) between the two, at the thalamus to achieve net sensory-motor adaptation. Thus, here we develop a system for the parallel combination of the input correlationbased learner (ICO) and the reward-based learner (actor-critic) as depicted in **Figure 1B**. The system works as a dual learner where the individual learning mechanisms run in parallel to guide the behavior of the agent. Both systems adapt their synaptic weights independently (as per their local synaptic modification rules) while receiving the same sensory feedback from the agent (environmental stimuli) in parallel. The final action that drives the agent is calculated as a weighted sum (**Figure 3** red circle) of the individual learning components. This can be described as follows:

$$
\rho\_{com}(t) = \xi\_{ico} o\_{ico}(t) + \xi\_{ac} o\_{ac}(t) \tag{1}
$$

where, *oico*(*t*) and *oac*(*t*) are the *t* time step outputs of the input correlation-based learner and the actor-critic reinforcement learner, respectively. *ocom*(*t*) represents the *t* time step combined action. The key parameters here that govern the learning behavior are the synaptic weights of the output neuron projection from the individual components (ξ*ico* and ξ*ac*). These govern the degree of influence of the two learning systems, on the net action of the agent. Previously, a simple and straight forward approach was undertaken in Manoonpong et al. (2013), where an equal contribution (ξ*ico* = ξ*ac* = 0.5) of ICO and actor-critic RL for controlling an agent was considered. Although this can lead to successful solutions in certain goal-directed problems, it is sub-optimal due to the lack of any adaptive balancing mechanism. Intuitively for associative learning problems with immediate rewards the ICO system learns quickly as compared to distal reward based goaldirected problems where, the ICO learner can provide guidance to the actor-critic learner. In particular depending on the type of problem, the right balance between the two learners needs to be achieved in an adaptive manner.

While there is evidence on the direct communication (Bostan et al., 2010) or combination of the subcortical loops from the cerebellum and the basal ganglia (Houk et al., 2007), a computational mechanism underlying this combination has not been presented, so far. Here we propose for the first time, an adaptive combination mechanism of the two components, modeled in the form of a reward modulated heterosynaptic plasticity (RMHP) rule, which learns the individual synaptic weights (ξ*ico* and ξ*ac*) for the projections from these two components. It is plausible that such a combination occurs at the VA-VL region of the motor thalamic nuclei which has both pallido-thalamic (basal ganglia) and cerebello-thalamic projections (Sakai et al., 2000). Furthermore, a few previous experimental studies (Desiraju and Purpura, 1969; Allen and Tsukahara, 1974) suggested that the individual neurons of the VL (nearly 20%) integrate signals from the basal ganglia and the cerebellum along with some weak cerebral inputs2. Based on biological evidence of dopaminergic projections at the thalamus from the basal ganglia circuit (García-Cabezas et al., 2007; Varela, 2014) as well as cerebellar projections to the thalamic ventro-latral nucleus (Bosch-Bouju et al., 2013) (see Figures 42–47 in Lisberger and Thach, 2013) we consider here that such dopaminergic projections act as the neuromodulatory signal and triggers the heterosynaptic plasticity (Ishikawa et al., 2013). A large number of such heterosynaptic plasticity mechanisms contribute toward a variety of neural processes involving associative learning and development of neural circuits in general (Bailey et al., 2000; Chistiakova and Volgushev, 2009). Although there is currently no direct experimental evidence of heterosynaptic plasticity at thalamic nuclei, it is highly plausible that such interactions could occur on synaptic afferents as observed in the amygdala and the hippocampus (Vitureira et al., 2012). Here, we use the instantaneous reward signal as the modulatory input in order to induce heterosynaptic changes at the thalamic junction. Similar approach have also been used in some previous theoretical models of reward modulated plasticity (Legenstein et al., 2008; Hoerzer et al., 2012). Although the dopaminergic projections from the VTA to the Mthal are primarily believed to encode a reward prediction error (RPE) signal (Schultz and Dickinson, 2000), there exists considerable diversity in the VTA neuron types with a subset of these dopaminergic neurons directly responding to rewards (Cohen et al., 2012). Similar variability has also been observed in the single DA neuron recordings from memory guided sacadic tasks performed with primates (Takikawa et al., 2004). This suggests that although most dopaminergic neurons respond to a reward predicting conditional simuli, some may not strictly follow the canonical RPE coding (Cohen et al., 2012). It is important to note that, within this model, it is equally possible to use the reward prediction error (TD error, Equation 12) and still learn the synaptic weights of the two components in a stable manner, however with a negligibly slower weight convergence due to continuous weight changes (see Supplementary Figure 1).

Based on this RMHP plasticity rule the ICO and actor-critic RL weights are learned at each time step as follows :

$$
\Delta \xi\_{ico}(t) = \eta r(t) [o\_{ico}(t) - \bar{o}\_{ico}(t)] o\_{ic}(t), \tag{2}
$$

$$
\Delta \xi\_{ac}(t) = \eta r(t)[o\_{ac}(t) - \bar{o}\_{ac}(t)]o\_{ico}(t). \tag{3}
$$

Here *r*(*t*) is the current time step reward signal received by the agent, while *o*¯*ico*(*t*) and *o*¯*ac*(*t*) denote the low-pass filtered version of the output from the ICO learner and the actor-critic learner, respectively. They are calculated as:

$$
\bar{o}\_{ico}(t) = 0.9\bar{o}\_{ico}(t-1) + 0.1o\_{ico}(t),
$$

$$
\bar{o}\_{ac}(t) = 0.9\bar{o}\_{ac}(t-1) + 0.1o\_{ac}(t). \tag{4}
$$

The plasticity model used here is based on the assumption that the net policy performance (agent's behavior) is influenced by a single global neuromodulatory signal. This relates to the dopaminergic projections to the ventra-lateral nucleus in the thalamus as well as connections from the amygdala which can carry reward related signals that influence over all action selection. The RMHP learning rule correlates three factors: (1) the reward signal, (2) the deviations of the ICO and actor-critic learner outputs from their mean values, and (3) the actual ICO and actor-critic outputs. The correlations are used to adjust their respective synaptic weights (ξ*ico* and ξ*ac*). Intuitively here the heterosynaptic plasticity rule can be also viewed as a homeostatic mechanism (Vitureira et al., 2012). Such that, the equation 2 tells the system to increase the ICO learners weights (ξ*ico*) when the ICO output is coincident with the positive reward, while the third factor (*oac*) tells the system to increase ξ*ico* more (or less) when the actor-critic learner weights (ξ*ac*) are large (or small), and vice versa for Equation 3. This ensures that overall the ratio of weight change of the two learning components occurs at largely the same rate. Additionally in order to prevent uncontrolled divergence in the learned weights, homeostatic synaptic normalization is carried out specifically as follows:

$$
\xi\_{ico}(t) = \frac{\xi\_{ico}(t)}{\xi\_{ico}(t) + \xi\_{ac}(t)},
$$

$$
\xi\_{ac}(t) = \frac{\xi\_{ac}(t)}{\xi\_{ico}(t) + \xi\_{ac}(t)}.\tag{5}
$$

This ensures that the synaptic weights always add up to one and 0 < ξ*ico*, ξ*ac* < 1. In general this plasticity rule occurs on a very slow time scale which is governed by the learning rate parameter η. Typically convergence and stabilization of weights are achieved by setting η much smaller compared to the learning rate of the two individual learning systems (ICO and actor-critic). To get a more detailed view of the implementation of the adaptive combinatorial learning mechanism, interested readers should refer to algorithm 2 in the Supplementary Material.

#### **2.2. INPUT CORRELATION MODEL OF CEREBELLAR LEARNING**

In order to model classical conditioning of adaptive motor reflexes3in the cerebellum, we use a model-free, correlation based, predictive control learning rule called input correlation learning (ICO) (Porr and Wörgötter, 2006). ICO learning provides a fast and stable mechanism in order to acquire and generate sensory predictions for adaptive responses based solely on the correlations between incoming stimuli. The ICO learning rule (**Figure 3** Right) takes the form of an unsupervised synaptic modification mechanism using the cross-correlation between the incoming predictive input stimuli (predictive here means that the signals occur early) and a single reflex signal (late occurring). As depicted in **Figure 3** right, cortical perceptual input in the form of predictive signals (CS) represents the mossy fiber projections to the cerebellum microcircuit, while the Climbing fiber projections from the inferior olive that modulates the synaptic weights in the

2It is also plausible that integration of activity arising in basal ganglia and cerebellum might take place in the thalamus nuclei other than the VL-VA, since pallidal as well as cerebellar fibers are known histologically to terminate not only in the VL-VA but also in other structures (Mehler, 1971).

3The reflex signal is typically a default response to an unwanted situation. This acts as the unconditional stimulus occurring later in time, than the predictive conditional stimulus.

deep cerebellar nucleus are depicted in a simplified form with the differential region (*d*/*dt*).

The goal of the ICO mechanism is to behave as a forward model system (Porr and Wörgötter, 2006) that uses the sensory CS to predict the occurrence of the innate reflex signal (external predefined feedback signaling unwanted scenarios), thus letting the agent to react in an anticipatory manner to avoid the basic reflex altogether. Based on a differential Hebbian learning rule (Kolodziejski et al., 2008) the synaptic weights in the ICO scheme are modified using heterosynaptic interactions of the incoming inputs, depending on their order of occurrence. In general, the plastic synapses of the predictive inputs get strengthened if they precede the reflex signal and are weakened if their order of occurrence is reversed. As a result, the ICO learning rule drives the behavior depending on the timing of correlated neural signals. This can be formally represented as,

$$
\rho\_{i\alpha}(t) = \rho\_0 \mathbf{x}\_0(t) + \sum\_{j=1}^{K} \rho\_j(t) \mathbf{x}\_j(t). \tag{6}
$$

Here, *oico* represents the output neuron activation of the ICO system driven by the superposition of the plastic K-dimensional predictive inputs *xj*(*t*) = *x*1(*t*), *x*2(*t*),..., *xK*(*t*) 4 (differentially modified) and the fixed innate reflex signal *x*0(*t*). The synaptic strength of the reflex signal is represented by ρ0 and is fixed to the constant value of 1.0 in order to signal innate response to the agent. Using the cross-correlations between the input signals, our differential Hebbian learning rule modifies synaptic connections as follows:

$$
\Delta \rho\_{\dot{\jmath}}(t) = \mu \mathbf{x}\_{\dot{\jmath}}(t) \frac{d}{dt} \mathbf{x}\_0(t). \tag{7}
$$

Here, μ defines the learning rate and is typically set to a small value to allow slow growth of synaptic weights with convergence occurring once the reflex signal *xo* = 0 (Porr and Wörgötter, 2006). Thus, ICO learning allows the agent to predict the primary reflex and successfully generate early, adaptive actions. However, no explicit feedback of goodness of behavior is provided to the agent and thus only an anticipatory response can be learned without the explicit notion of how well the action allows reaching a desired (rewarding) goal location. As depicted in **Figure 3**, the output from the ICO learner is directly fed into the RMHP unit envisioned to be part of the ventro-lateral thalamic nucleus (Akkal et al., 2007; Bosch-Bouju et al., 2013).

#### **2.3. ACTOR-CRITIC RESERVOIR MODEL OF BASAL-GANGLIA LEARNING**

TD learning (Sutton, 1988; Suri and Schultz, 2001), in the framework of actor-critic reinforcement learning (Joel et al., 2002; Wörgötter and Porr, 2005), is the most established computational model of the basal ganglia. As explained in the previous section, the TD learning technique is particularly well suited for replicating or understanding how reward related information is formed and transferred by the mid-brain dopaminergic activity.

The model consists of two sub-networks, namely, the adaptive critic (**Figure 3** left, bottom) and the actor (**Figure 3** left, above). The critic is adaptive in the sense that it learns to predict the weighted sum of future rewards taking into account the current incoming sensory stimuli and the actions (behaviors) performed by the agent within a particular environment. The difference between the predicted "value" of sum of future rewards and the actual measure acts as the temporal difference (TD) prediction error signal that provides an evaluative feedback (or reinforcement signal) to drive the actor. Eventually the actor learns to perform the proper set of actions (policy5 ) that maximize the weighted sum of future rewards as computed by the critic. The evaluative feedback (TD error signal) in general acts as a measure of goodness of behavior that, overtime, lets the agent learn to anticipate reinforcing events. Within this computational framework, the TD prediction error signal and learning at the critic are analogous to the dopaminergic (DA) activity and the DA dependent long term synaptic plasticity in the striatum (**Figure 2B**), while the remaining parts of striatal circuitry can be envisioned as the actor which uses the TD modulated activity to generate actions, which drives the agent's behavior.

Inspired by the reservoir computing framework (Maass et al., 2002; Jaeger and Haas, 2004), here we use a chaotic random recurrent neural network (RNN) (Sussillo and Abbott, 2009; Rajan et al., 2010) as the adaptive critic (cortico-striatal circuitry and the DA system) connected to a feed-forward neural network, serving the purpose of the part of striatum that performs action selection (Gurney et al., 2001) and then relays it to the motor thalamus via projections from the globus pallidus and substantia nigra. This provides an effective framework to model a continuous actorcritic reinforcement learning scheme, which is particularly suited for goal-directed learning in continuous state-action problems, while at the same time maintaining a reasonable level of biological abstraction (Fremaux et al., 2013). Here, the reservoir network can be envisioned as analogous to the cortex and its inherent recurrent connectivity structure, and the readout neurons serving as the striatum, with plastic projections from the recurrent layer, as the modifiable cortico-striatal connections (Hinaut and Dominey, 2013). The reservoir network is constructed as a generic network model of *N* recurrently connected neurons with high sparsity (refer to Supplementary Material for details) and fixed synaptic connectivity. The connections within the recurrent layer are drawn randomly in order to generate a sparsely connected network of inhibitory and excitatory synapses. A subset of the reservoir neurons receive input connections (fixed synaptic strengths) as external driving signals and has an additional output layer of neurons that learns to produce a desired response based on synaptic modification of weights from the reservoir to output neurons. The input connections along with the large recurrently connected reservoir network represents the main cortical microcircuit-to-striatum connections, while the output layer neural activity can be envisioned as striatal neuronal responses. In this case, the reservoir critic provides an input (sensory stimuli) driven dynamic network with a large repertoire of signals

4This x(t) is different from the neural state activation vector **x**(*t*) of Equation 9.

5In reinforcement learning, policy refers to the set of actions performed by an agent that maximizes it's average future reward.

that is used to predict the value function *v* (average sum of future rewards). *v*(*t*) approximates the accumulated sum of the future rewards *r*(*t*) with a given discount factor γ (0 ≤ γ < 1)6 as follows:

$$\nu(t) = \sum\_{i=1}^{\infty} \nu^{i-1} r(t+i). \tag{8}$$

In our model, the membrane potential at the soma (at time *t*) of the reservoir neurons, resulting from the incoming excitatory and inhibitory synaptic inputs, is given by the *N* dimensional vector of neuron state activation's, **x**(*t*) = *x*1(*t*), *x*2(*t*),..., *xN*(*t*). The input to the reservoir network, consisting of the agent's states (sensory input stimuli from the cerebral cortex), is represented by the *K* dimensional vector **u**(*t*) = *u*1(*t*), *u*2(*t*),..., *uK*(*t*). The recurrent neural activity within the dynamic reservoir varies as a function of its previous state activation and the current driving input stimuli. The recurrent network dynamics is given by,

$$\tau \dot{\mathbf{x}}(t) = -\mathbf{x}(t) + \mathbf{g} \mathbf{W}\_{\text{sys}} \mathbf{z}(t) + \mathbf{W}\_{\text{in}} \mathbf{u}(t) + \mathbf{b},\tag{9}$$

$$
\hat{\nu}(t) = \tanh(\mathbf{W}\_{out}\mathbf{z}(t)),
\tag{10}
$$

$$z\_i(t) = \tanh(\alpha x\_i(t) + \beta). \tag{11}$$

The parameters **W***in* and **W***sys* denote the input to reservoir synaptic weights and the recurrent connection weights within the reservoir, respectively. The parameter *g* (Sompolinsky et al., 1988) acts as the scaling factor for the recurrent connection weights allowing different dynamic regimes from stable to chaotic being present in the reservoir. Similar to Sussillo and Abbott (2009) we select *g* such that the network exhibits chaotic dynamics as spontaneous behavior before learning and maintains stable dynamics after learning, with the help of feedback connections and neuronal activation homeostasis via intrinsic plasticity (Triesch, 2005; Dasgupta et al., 2013a). The RNN does not explicitly model action potentials, but describes neuronal firing rates, where in, the continuous variable *zi* is the instantaneous firing rate of the reservoir neurons and is calculated as a non-linear saturating function of the state activation *xi* (Equation 11). The output layer consists of a single neuron whose firing rate *v*ˆ(*t*) is calculated at time *t* based on equation 10, as a non-linear transformation of the weighted projections of the internal reservoir neuron firing rates **z**(*t*). Here the parameter **W***out* denotes the *N* × *K* dimensional reservoir to output connection synaptic weights. Each unit in the network also receives a constant bias signal *bi*, represented in equation 9 as the *N* dimensional vector **b**. The overall time scale of the RNN and the leak rates of individual reservoir neurons are controlled by the parameter τ .

Based on the TD learning principle, the primary goal of the reservoir critic is to predict *v*(*t*) such that the TD error δ is minimized over time. At each time point *t*, δ is computed from the current (*v*ˆ(*t*)) and previous (*v*ˆ(*t* − 1)) value function predictions (reservoir output), and the current reward signal *r*(*t*), as follows:

$$
\delta(t) = r(t) + \mathcal{y}\hat{\boldsymbol{\nu}}(t) - \hat{\boldsymbol{\nu}}(t-1). \tag{12}
$$

The output weights **W***out* are calculated using the recursive least squares (RLS) algorithm (Haykin, 2002) at each time step, while the sensory stimuli **u**(*t*) are being fed into the reservoir. **W***out* are calculated such that the overall TD-error (δ) is minimized. We implement the online RLS algorithm using a fixed forgetting factor (λ*RLS* < 1) as given in **Algorithm 1**.

As proposed in Triesch (2005) and Dasgupta et al. (2013a) we implement a generic intrinsic plasticity mechanism based on the Weibull distribution for unsupervised adaptation of the reservoir neuron non-linearity using a stochastic decent algorithm to adapt the scale α and shape parameters β of the saturating function in Equation 11. This allows the reservoir to homoeostatically maintain a stable firing rate while at the same time it drives the neuron activities to a non-chaotic regime. It is also important to note that one of the primary assumptions of the basic TD learning rule is a Markovian one, which considers future sensory cues and rewards depending only on the current sensory cue without any memory component. The use of a reservoir critic (due to the inherent fading temporal memory of recurrent networks Lazar et al., 2007) breaks this assumption. As a result, such design principle extends our model to problems with short term dependence of immediate sensory stimuli on the preceding history of stimuli and reward (see **Figure 4** for a simulated example of local temporal memory in reservoir neurons).

The actor (**Figure 3** left above) is designed as a single stochastic neuron, such that for a one dimensional action generation the output (*Oac*) is given as:

#### **Algorithm 1: Online RLS algorithm for learning reservoir to output neuron weights.**

*Initialize*: **W***out* = 0, exponential forgetting factor (λ*RLS*) is set to a value less than 1 (we use 0.85) and the auto-correlation matrix ρ is initialized as ρ(0) = **I**/β, where **I** is unit matrix and β is a small constant.

*Repeat*: At time step *t*

Step 1: For each input signal **u**(*t*), the reservoir neural firing rate vector **z**(*t*) and network output *v*ˆ(*t*) are calculated using equation 11 and equation 10.

Step 2: Online error *e*(*t*) calculated as: *e*(*t*) ← δ(*t*)

Step 3: Gain vector **K**(*t*) is updated as: **K**(*t*) ← ρ(*t*−1)**z**(*t*) λ*RLS*+**z***T* (*t*)ρ(*t*−1)**z**(*t*)

Step 4: Update the auto-correlation matrix ρ(*t*) ρ(*t*) ← 1 λ*RLS* ρ(*t* − 1) − *K*(*t*)**z***T*(*t*)ρ(*t* − 1)

Step 5: Update the instantaneous output weights **W***out*(*t*) **W***out*(*t*) ← **W***out*(*t* − 1) + *K*(*t*)*e*(*t*)

Step 6: *t* ← *t* + 1

*Until*: The maximum number of time steps is reached.

6The discount factor helps assigning decreasing value to rewards further away in the past as compared to the current reward.

*oac*(*t*) = (*t*) +- *K wi*(*t*)*ui*(*t*), (13)

(not shown here). Spontaneous dynamics then unfolds in the system based on Equation 9. The lower right panel plots the activity of 5

*i* = 1

where *K* denotes the dimension (total number) of sensory stimuli (**u**(*t*)) to the agent being controlled. The parameter *wi* denotes the synaptic weights for the different sensory inputs projecting to the actor neuron. Stochastic noise is added to the actor via (*t*), which is the exploration quantity updated at every time step. This acts as a noise term, such that initially exploration is high, and the agent needs to navigate the environment more if the expected cumulative future reward *v*(*t*) is sub-optimal. However, as the agent learns to successfully predict the maximum cumulative reward (value function) over time, and the net exploration is decreased. As a result (*t*) gradually tends toward zero as the agent starts to learn the desired behavior (correct policy). Using Gaussian white noise σ (zero mean and standard deviation one) bounded by the minimum and maximum limits of the value function (*vmin* and *vmax*), the exploration term is modulated as follows:

$$\epsilon(t) = \Omega \sigma(t) \cdot \min\left[0.5, \max\left(0, \frac{\nu\_{\max} - \hat{\nu}(t)}{\nu\_{\max} - \nu\_{\min}}\right)\right]. \tag{14}$$

Here, is a constant scale factor selected empirically (see Supplementary Material for details). The actor learns to produce the correct policy, by an online adaptation (**Figure 3** left above) of its synaptic weights *wi* at each time step as follows:

$$
\Delta \boldsymbol{w}\_i(t) = \boldsymbol{\pi}\_a \boldsymbol{\delta}(t) \boldsymbol{u}\_i(t) \boldsymbol{\epsilon}(t), \tag{15}
$$

where τ*a* is the learning rate such that 0 < τ*a* < 1. Instead of using direct reward *r*(*t*) to update the input to actor neuron synaptic weights, using the TD-error (i.e., error of an internal reward) allows the agent to learn successful behavior, even in cases of delayed reward scenarios (reward is not given uniformly for each time step but is delivered as a constant value after a set of actions were performed to reach a specific goal). In general, once the agent learns the correct behavior, the exploration term ((*t*)) should become zero, as a result of which no further weight change (Equation 15) occurs and *oac*(*t*) represents the desired action policy, without any additional noise component.

the network exhibits considerable fading memory of the brief

#### **3. RESULTS**

incoming input stimuli.

In order to test the performance of our bio-inspired adaptive combinatorial learning mechanism, and validate the interaction through sensory feedback, between reward-based learning (basal ganglia) and correlation-based learning (cerebellum) systems, we employ a simulated, goal-directed decision making scenario of foraging behavior. This is carried out within a simplified paradigm of a four-wheeled robot navigating an enclosed environment, with gradually increasing task complexity.

#### **3.1. ROBOT MODEL**

The simulated wheeled robot NIMM4 (**Figure 5**) consists of a simple body design with four wheels whose collective degree of rotation controls the steering and the over all direction of motion. It is provided with two front infrared sensors (*IR*1 and *IR*2) which can be used to detect obstacles to its left or right side, respectively. Two relative orientation sensors (μ*G* and μ*B*) are also provided, which can continuously measure the angle of deviation of the robot with respect to the green (positive) and blue (negative) food sources. They are calibrated to take values in the interval [−180◦, 180◦] with the angle of deviation μ*G*,*B* = 0*o* when the respective goal is directly in front of the robot, μ*G*,*B* is positive when the goal locations are to the right of the robot and negative for the opposite case. In addition NIMM4 also consists of two relative position sensors (*DG*,*B*) that can calculate it's relative straight line distance to a goal, taking values in the interval [0, 1], with the

**FIGURE 5 | Simulated mobile robot system for goal-directed behavior task. (Top)** The mobile robot NIMM4 with different types of sensors. The relative orientation sensor μ is used as state information for the robot. **(Bottom)** Variation of the relative orientation μ*G* to the green goal. the front left and right infrared sensors *IR*1 and *IR*2 are used to detect obstacles in front of the robot. Direction control for the robot is maintained using the quantity *Usteering* calculated by the individual learning components (ICO and actor-critic) and then fed to the robot wheels to generate forward motion or steering behavior. Sensors *DG* and *DB* measure straight line distance to the goal locations.

respective sensor reading tending to zero, as the robot gets closer to the goal location and vice versa.

#### **3.2. EXPERIMENTAL SETUP**

The experimental setup (**Figure 6**) consists of a bounded environment with two different food sources (desired vs punishing) located at fixed positions. The primary task of the robot is to navigate the environment such that, eventually, it should learn to steer toward the food source that leads to positive reinforcements (green spherical ball in **Figures 6A–C**) while avoiding the goal location that provides negative reinforcements or punishments (blue spherical ball), within a specific time interval. The main task is designed as a continuous state-action problem with a distal reward setup (Reinforcement zone in **Figure 6**), such that the robot starts at a fixed spatial location with random initial orientation ([−60◦, 60◦]) and receives the positive or negative reinforcement signal only within a radius of specific distance (*DG*,*B* = 0.2) from the two goal locations. Within this boundary, for the green goal it receives a continuous reward of +1 at every time step and a continuous punishment of −1 in case of the blue goal, respectively. At other locations along the environment no reinforcement signal is given to the robot.

The experiments are further divided into three different scenarios of, foraging without an obstacle (case I), with single obstacle (case II) and a dynamic foraging scenario (case III), demonstrating different degrees of reward modulated adaptation between the two learning systems in different environments. In all scenarios, the robot can continuously sense its angle of deviation to the two goals with μ*G*,*B* always active. This acts as a Markov decision process (MDP) such that, the next sensory state of the robot depends on the sensory information for the current state of the robot and the action it performed, and is conditionally independent of all the previous sensory states and actions. Detecting the obstacle results in negative reinforcement (continuous −1 signal) triggered by the front infrared sensors (*IR*1,2 > 1.0). Furthermore, hitting the boundary wall in the arena results

#### **FIGURE 6 | Three different scenarios for the goal-directed foraging**

**task. (A)** Environmental setup without an obstacle case. Green and Blue objects represent the two food sources with positive and negative rewards, respectively. The red dotted circle indicates the region where the turning reflex response (from the ICO learner) kicks in. The robot is started from and reset to the same position, with random orientation at the beginning of each trial episode. **(B)** Environmental setup with an obstacle. In addition to the previous setup, a large obstacle is place in

the middle of the environment. The robot needs to learn to successfully avoid it and reach the rewarding food source. Collisions with the obstacle (triggered by *IR*1 and *IR*2) generate negative rewards (−1 signal) to the robot. **(C)** Environmental setup with dynamic switching of the two objects. It is an extended version of the first scenario. After every 50 trials the reward zones are switched such that the robot has to dynamically adjust to the new positively reinforced location (food) and learn a new trajectory from the starting location.

in a negative reinforcement signal (−1), with the robot being reset to the original starting location. Although the robot is provided with relative distance sensors, sensory stimuli (state information) is provided using only the angle of deviation sensors and the infrared sensors. The reinforcement zone (distance of *DG*,*B* = 0.2) is also used as the zone of reflex to trigger a reflex signal for the ICO learner. Fifty runs were carried out for each setup in all cases. Each run consisted of a maximum of 150 trials. The robot was reset if the maximum simulation time of 15 s was reached, or if it reaches one of the goal locations or if it hits a boundary wall, which ever occurs earlier.

#### **3.3. CEREBELLAR SYSTEM: ICO LEARNING SETUP**

The cerebellar system in the form of ICO learning (**Figure 3** right) was setup as follows: μ*G*,*B* were used as predictive signals (CS). Two independent reflex signals (*x*0,*B* and *x*0,*G*, see equation 6) were configured with one for blue food source and the other for the green food source (US). The setup was designed following the principles of delayed conditioning experiments, where, an overlap between the CS and the US stimuli needs to exist in order for the learning to take place. The reflex signal was designed (measured in terms of the relative orientation sensors of the robot) to elicit a turn toward a specific goal once the robot comes within the reflex zone (inside the dotted circle in **Figures 6B,C**). Irrespective of the kind of goal (desired or undesired) the reflex signal drives the robot toward it with a turn proportional to the deviations defined by μ*G*,*B* i.e., large deviations cause sharper turns. The green and the blue ball were placed such that there was no overlap between the reflex areas, hence only one reflex signal per goal, got triggered at a time. In other words, the goal of the ICO learner is simply to learn to steer toward a food location without any knowledge of it's worth. This is representative of an adaptive reflexive behavior as observed in rodent foraging studies where in the behavior is guided without explicit rewards, but just driven by conditioning between the CS-US stimuli, such that the robot or animal learns to favor certain spots in the environments without any knowledge of their worth. The weights of the ICO learner ρμ*G* and ρμ*B* (Equation 6) with respect to the green and blue goals were initialized to 0.0. If the positive derivative of the reflex signal becomes greater than a predefined threshold, the weights change and otherwise they remain static, i.e., a higher change in ρμ*G* in comparison to ρμ*B* would mean that the robot gets drawn toward the green goal more.

#### **3.4. BASAL GANGLIA SYSTEM: RESERVOIR ACTOR-CRITIC SETUP**

The basal ganglia system in the form of a reservoir based actorcritic learner was setup such that, the inputs to the critic and actor networks (**Figure 3** left) consisted of the two relative orientation sensor data μ*G* and μ*B* and the front left and right infrared sensors (*IR*1 and *IR*2) of the robot (**Figure 4**). Although the robot also contains relative distance sensors, these were not used as state information inputs. This makes the task less trivial, such that sufficient but not complete information was provided to the actor-critic RL network. The reservoir network for the critic consisted of *N* = 100 neurons and one output neuron that estimates the value function *v*(*t*) (Equation 10). Reservoir input weights *Win* were drawn from an uniform distribution [−0.5, 0.5] while the reservoir recurrent weights *Wsys* were drawn from a Gaussian distribution of mean 0 and standard deviation *g*2/*N* (see Equation 9). Here *g* acts as the scaling factor for *Wsys*, and it was designed such that there is only 10% internal connectivity in *Wsys* with a scaling factor of 1.2. The reward signal *r*(*t*) (Equation 12) was set to +1 when the robot comes close (reflex/reinforcement zone) to the green ball and to −1 when it comes close to the blue ball. A negative reward of −1 was also given for any collisions with the boundary walls or obstacle. At all other locations within the environment, the robot receives no explicit reward signal. Thus, the setup is designed keeping a delayed reward scenario in mind, such that earlier actions lead to a positive or negative reward, only when the robot enters the respective reinforcement/reflex zone. The synaptic weights of the actor with respect to the two orientation sensors (*w*μ*G* and *w*μ*B* ) were initialized to 0.0, while the weights with respect to the infrared sensors (*wIR*1 and *wIR*2 ) were initialized to 0.5 (equation 13). After learning, a high value of *w*μ*G* and a low value of *w*μ*B* would drive the robot toward the green goal location and away from the blue goal. The weights of the infrared sensor inputs effectively control the turning behavior of the robot when encountered with an obstacle (higher *wIR*1—right turn, higher *wIR*2—left turn). The parameters of the adaptive combinatorial network are summarized in the Supplementary Tables 1–3.

#### **3.5. CASE I: FORAGING WITHOUT OBSTACLE**

In the simplest foraging scenario the robot was placed in an environment with two possible food sources (green and blue) and without any obstacle in between (**Figure 6A**). In this case the green food source provided positive reward while the blue food source provided negative reward. The goal of the combined learning mechanism was to make the robot successfully steer toward the desired food source. **Figure 7A** shows simulation snapshots of the behavior of the robot as it explores the environment. As observed from the trajectory of the robot, initially it performed a lot of exploratory behavior and randomly moved around in the environment, but eventually it learned to move solely toward the green goal. This can be further analyzed looking at the development of the synaptic weights of the different learning components as depicted in **Figure 8**. As observed in **Figure 8C** due to the simple correlation mechanism of the ICO learner (cerebellar system), the ICO weights adapt relatively faster as compared to the actor. Due to random explorations (**Figure 9B**) in the beginning, in the event of the blue goal being visited more frequently, reflexive pull toward blue goal - ρμ*B* is greater than toward the green goal ρμ*G* . However, after sufficient explorations, as the robot starts reaching the green goal more frequently, ρμ*G* also starts developing. This is counteracted by the actor weights (basal ganglia system), where in, there is a higher increase in *w*μ*G* (orientation sensor input representing angle of deviation from green goal) as compared to *w*μ*B* (orientation sensor input representing angle of deviation from blue goal). This is caused as result of the increased positive rewards received from the green goal (**Figure 9A**) that causes the TD-error to modulate the actor weights (equation 15) accordingly. At the same time no significant change is seen in the infrared sensor input weights (**Figure 8B**), due to the fact that in this scenario, the infrared sensors get triggered only on collisions

with the boundary wall and remain dormant otherwise. Recall that the infrared sensor weights were initialized to 0.5.

Over time as the robot moves more toward the desired food source, the ICO weights also stabilize with the reflex toward the green goal being much stronger. This also leads to a reduction of the exploration noise (**Figure 9B**), and the actor weights eventually converge to a stable value (**Figures 8A,B**). Here, the slow RMHP rule performs a balancing act between the two learning systems with initial higher weight of the actor-critic learner and then a switch toward the ICO system, once the individual learning rules have converged. **Figure 9C** shows the development of the value function (*v*(*t*)) at each trial, as estimated by the critic. As observed initially the critic underestimates the total value due to high explorations and random navigation in the environment. However, as the different learning rules converge, the value function starts to reflect the total accumulated reward with stabilization after 25 trials (each trials consisted of approximately 1000 time steps).

This is also clearly observed from the change of the orientation sensor readings shown in **Figure 9D**. Although there is considerable change in the sensor readings initially, after learning, the orientation sensor toward the green goal (μ*G*) records positive angle, while the orientation from the blue goal μ*B* records considerably lower negative angles. This indicates that the robot learns to move stably toward the positively rewarded food source and away from the oppositely rewarded blue food source. Although this is the simplest foraging scenario, the development of the RMHP weights ξ*ico* and ξ*ac* (**Figure 8D**) depicts the adaptive combination of the basal gangliar and cerebellar learning systems for goal-directed behavior control. Here the cerebellar system (namely ICO) acts as a fast adaptive reflex learner that guides and shapes the behavior of the

**FIGURE 8 | Synaptic weight change curves for the static foraging tasks without obstacle and with single obstacle. (A)** Change in the synaptic weights for actor-critic RL learner. Here *w*μ*G* corresponds to the input weights of the orientation sensor toward the green goal and *w*μ*B* corresponds to the input weights of the orientation sensor toward the blue goal. **(B)** Change in the weights of the two infrared sensor inputs of the actor. *wIR*1 is the left IR sensor weight, *wIR*2 is the right IR sensor weights. **(C)** Change in the synaptic weights of the ICO learner. ρμ*G* is the CS stimulus weight for the orientation sensor toward green, ρμ*B* the CS stimulus weight for the orientation sensor

toward blue. **(D)** Learning curve of the RMHP combined learning mechanism showing the change in the weights of the ICO network output (depicted in red). ξ*ico* is weight of the ICO network output. ξ*ac* is weight of the actor-critic RL network output (depicted in black). **(E–H)** Show the change in the weights corresponding to the single obstacle static foraging task. In all the plots the gray shaded region marks the region of convergence for the respective synaptic weights. Three different timescales exist in the system, with the ICO learning being the fastest, actor-critic RL being intermediate and the adaptive combined learning being the slowest. (see text for more details.)

Dasgupta et al. Neuromodulated adaptive neural combinatorial learning

reward-based learning system. Although both the individual systems eventually converge to provide the correct weights toward the green goal, the higher strength of the ICO component (ξ*ico*) leads to a good trajectory irrespective of the starting orientation of the robot. This is further illustrated in the simulation video showing three different scenarios of only ICO, only actor-critic and the combined learning cases, see Supplementary Movie 1.

#### **3.6. CASE II: FORAGING WITH SINGLE OBSTACLE**

In order to evaluate the efficacy of the two learning systems and their cooperative behavior, the robot was now placed in a slightly modified environment (**Figure 6B**). As in the previous case, the robot still starts from a fixed location with initial random orientations. However, it now has to overcome an obstacle placed directly in front (field of view), in order to reach the rewarding food source (green goal). Collisions with the obstacle, during learning, resulted in negative rewards (−1) triggered by the front left (*IR*1) and right (*IR*2) infrared sensors. This influenced the actor-critic learner to modulate the actor weights via

exploration steadily decreases and the value function prediction also reaches near convergence at 25 trials (1 trial approximates 1000 time steps). The thick black line represents the average value calculated over 50 runs of the experiment with standard deviation given by the shaded region. **(D)** Plots of the two orientation sensor readings (in degrees) for the green (μ*G*) and the blue (μ*B*) goals, averaged over 50 runs. During initial exploration the angle of the deviation of the robot from the two goals changes randomly. However, after convergence of the learning rules, the orientation sensor readings stabilize with small positive angle of deviation toward the green goal and large negative deviation from the blue goal. This shows that post learning, the robot steers more toward the green goal and away from the blue goal. Here the thick lines represent average values and the shaded regions represent standard deviation.

TD-error and generate turning behavior around the obstacles. In parallel, the ICO system, still learns only a default reflexive behavior of getting attracted toward either of the food sources by a magnitude proportional to its proximity to them (same as case I), irrespective of the associated rewards. As observed from the simulation snapshots in **Figure 7B**, after initial random exploration, the robot learns the correct trajectory to navigate around the obstacle and reach the green goal. From the synaptic weight development curves for the actor neuron (**Figure 8E**) it is clearly observed that although initially there is a competition between *w*μ*G* and *w*μ*B* , after sufficient exploration, as the robot gets more positive rewards by moving to the green food source, the *w*μ*G* weight becomes larger in magnitude and eventually stabilizes.

Concurrently in **Figure 8F**, it can be observed that unlike the previous case the left infrared sensor input weight *wIR*1 gets considerably higher as compared to *wIR*2 . This is indicative of the robot learning the correct behavior of turning right in order to avoid the obstacle and reach the green goal. However, interestingly, as opposed to the simple case (no obstacle) the ICO learner tries to pull the robot more toward the blue goal, as seen from the weight development of ρμ*G* and ρμ*B* in **Figure 8G**. This behavior can be attributed to the fact that, as the robot reaches the blue object in the beginning, the fast ICO learner provides high weights for a reflexive pull toward the blue as opposed to the green goal. As learning proceeds and the robot learns to move toward the desired location (driven by the actor-critic system), the ρμ*G* weight also increases, however it still continues to favor the blue goal. As a result in order to learn the correct behavior the combined learning systems needs to favor the actorcritic mechanism more as compared to the naive reflexives from the ICO. This is clearly observed from the balancing between the two as depicted in the ξ*ico* and ξ*ac* weights in **Figure 8H**. Following the stabilization of the individual learning system weights, the combined learner provides much higher weighting of the actor-critic RL system. Thus, in this scenario, due to the added complexity of an obstacle, one sees that the reward modulated plasticity (RMHP rule) learns to balance the two interacting learning systems, such that the robot still performs the correct decisions overtime (see the simulation run from Supplementary Movie 2).

## **3.7. CASE III: DYNAMIC FORAGING (REVERSAL LEARNING)**

A number of modeling as well as experimental studies of decision making (Sugrue et al., 2004) have considered the behavioral effects of associative learning mechanisms on dynamic foraging tasks as compared to static ones. Thus, in order to test the robustness of our learning model, we changed the original setup (**Figure 6C**), such that, initially a positive reward (+1) is given for the green object and a negative reward (−1) for the blue one. This enables the robot to learn moving toward the green object while avoiding the blue object. However, after every 50 trials the sign of the rewards was switched such that now the blue object received positive reward, and the green goal the opposite. As a result the learning system needs to quickly adapt to the new situation and learn to navigate to the correct target. As observed in the **Figure 10B** initially the robot performs random explorations receiving a mixture of positive and negative rewards, however after sufficient trials, the robot reaches a stable configuration (exploration drops to zero) and receives positive rewards concurrently (**Figure 10A**). This corresponds to the previous case of learning to move toward the green goal. As the rewards were switched, the robot then obtained negative reward when it moved to the green object. As a consequence, the exploration gradually increased again; thereby the robot also exhibited random movements. After successive trials, a new stable configuration was reached with the exploration dropping to zero and now the robot received more positive rewards, however for the other target (blue object). This is depicted with more clarity, in the simulation snapshots in **Figure 7C** (beginning—random explorations, learn 1—reaching green goal, learn 2—reaching blue goal).

In order to understand how the combined learning mechanism handles this dynamic switching, in **Figure 11** we plot the synaptic weight developments of the different components.

Initially the robot behavior is shaped by the ICO weights (**Figure 11B**) which learn to steer the robot to the desired location, such that the reflex toward green object (ρμ*G* ) is stronger than that toward the blue object (ρμ*B* ). Furthermore, as the robot