The CRL is based on the assumption that at every synapse, two synaptic eligibility traces, one for LTP and one for LTD, can simultaneously be activated through the temporal firing patterns of the pre- and post-synaptic cells in a Hebbian manner. Moreover, the activity of these traces can saturate at different levels and in the absence of any cellular activity their activation level decays with different time constants. These assumptions of the temporal dynamics of these traces can be formulated mathematically as a pair of coupled differential equations of the form:

where the labels p and d refer to LTP and LTD synaptic eligibility traces, respectively. Thus, for a ∈ {p, d}, Tija represents the activation of the synaptic eligibility trace at a synapse labeled by the index ij, (i.e., between the j-th presynaptic and i-th postsynaptic cells), τ_a is the characteristic decay time constant, and Tmaxa sets the saturation level for the activation of the trace. H^a stands for the Hebbian rule, which itself depends on the firing rates or spike patterns of the pre- and post-synaptic cells S_j(t) and S_i(t), respectively. In general there can be different Hebbian rules for every synapse and thus we can have Ha(Si,Sj)=Hija. In many applications we will use a rate-dependent Hebbian rule of the form Hija=SiSj, where the S_k, the average firing rates, are computed using Equation (25). The above equations can be cast in a more familiar form

where the indices i, j have been dropped, τ~a is the effective time constant and T~a is the effective saturation activity of the trace given by

The general solution of Equation (3) for a time-dependent Hebbian term, H^a = H^a(t), is

where T^a(0) is the activity of the trace at time t = 0.

In order to gain some insight into the dynamics of the traces we consider the special case of a constant Hebbian term, i.e., H^a = const. Assuming that H^a is non-zero during a period of time (rising phase) and then turns off (falling phase), the dynamics of the traces is as follows. When H^a is non-zero, the trace dynamics will converge exponentially to the non-zero level.

Thus, during the rising phase (H^a > 0), the dynamics of the eligibility trace take the form

where R labels this solution for the rising phase.

During the falling phase (H^a = 0), the trace will decay to zero with a time constant τ^a following an exponential decay (Figure 1A)

where F labels this solution for the falling phase and TRa(tstim) is the activity of the trace at the end of the rising phase (t = t_stim). The time course of the trace activity for this example is illustrated in Figure 1A.

Note that these equations are based on the simplifying assumption that H^a is constant throughout the different phases. However, in general the Hebbian term H^a will vary in time and the dynamics of the traces will deviate from the simple description given earlier. Since in truth H^a will fluctuate over time and from trial to trial, the equations describing the activity of the traces characterize the average behavior.

An example of the temporal dynamics of the activity of the trace for a smoothly time-varying Hebbian term is given in Figure 1B. In this example, the activity of the trace increases rapidly in a non-exponential way during the rising phase. As the magnitude of the Hebbian term starts decreasing, the initial decay of the trace is slow but as the magnitude of the Hebbian term decreases to zero (gray shaded area), the activity follows an exponential decay, illustrated by the green-dotted line, with its characteristic decay constant τ^a. Further examples from actual simulation are given in Figure 2C.

Throughout this paper we use the simplified expression Hija=ηaHij, to indicate that the Hebbian terms for LTP and LTD are identical up to a multiplicative constant. Under this assumption Equations (4, 5) are transformed to:

with time constant

Given these synaptic eligibility traces, we will further assume that synaptic efficacies are only modified in the presence of time dependent reward signals, which are most likely to be neuromodulators. This assumption can be simply expressed by the equation:

where L_ij is the synaptic efficacy between neurons i and j, and R^p, R^d are the time dependent reward signals that convert the eligibility traces for LTP and LTD, respectively, to synaptic efficacy and η is the learning rate. For many properties of the learning rules a single neuromodulator is sufficient, in which case R^p = R^d. The temporal profile of the neuromodulator release is assumed to be shorter than the duration of the trial, T_trial.

The solution of this equation has the form:

where L_ij(0) is the synaptic efficacy at the beginning of the trial. Steady-state is obtained when the synaptic efficacy no longer changes, i.e., L_ij(T_trial) = L_ij(0).

We will consider several simple cases for the temporal profile of neuromodulators. The simplest model, which was used previously (Huertas et al., 2015), assumes that Rp(t)=Rd(t)=Rδ(t-treward) where t_reward is the time that the reward signal arrives in the cortical area. When examining the effects of different neuromodulators, we will use a model of the form Rp(t)=R0pδ(t-treward) and Rd(t)=R0dδ(t-treward) where R0d≠R0p. Another possibility considered here is that the neuromodulator level is constant within an time interval between t₁ and t₂. Obviously these are very simple equations, but they can be replaced with more biophysically realistic ones without a qualitative change in the behavior of these rules. In general the two neuromodulators might have different release profiles and might not only depend of reward, but possibly also on expected reward or on action.

As described previously by He et al. (2015) and as will be described in the following sections, we find that in order to obtain stable reinforcement learning, the saturation level of the LTD associated trace must be larger than the saturation of the LTP trace (Tmaxd>Tmaxp). Further, to obtain learning with a delayed reward in a recurrent network, the decay of the LTP trace must be slower than that of the LTD trace (τ_p > τ_d). We also show below that to obtain stable learning in a feedforward network, the rise time of the LTP trace must be faster than the rise time of the LTD trace (τ~p<τ~d).

We characterize the consequences of CRL in a network of integrate-and-fire neurons, as described in the Methods section, using the parameters in Table 1. We have previously shown (He et al., 2015) that with this rule it is possible to train the network to report reward times, in slowly decaying sustained neurons in V1 (Shuler and Bear, 2006). In Figure 2 we show a summary of these results. It shows that the network is able to learn spike timing for a broad range of reward times in a noise free network composed of only excitatory neurons (Figure 2A), as shown previously (He et al., 2015). We also show similar results in a network composed of both inhibitory and excitatory neurons and with realistic levels of variability (Figure 2B).

The process of training a recurrent network with the two-trace model is shown Figure 2C. Initially (Figure 2C, top left) the network activity decays quickly to baseline after the end of the stimulus, and at the time of reward the LTP related trace (Figure 2 bottom left) is larger than the LTD trace, due to the longer time constant of the LTP trace. This results in a potentiation of the synaptic efficacies, leading to longer lasting network activity (Figure 2C top, center). During learning, at the time of reward, the LTP trace is still larger than the LTD trace (Figure 2C bottom, center), leading to further enhancement of synaptic efficacies. Finally when the network activity approaches the time of reward (Figure 2C top, right) the LTP and LTD traces become equal (Figure 2C bottom, right), the synaptic efficacies no longer change, and the state of the network is stabilized.

Due to the strong connection between neurons, and in particular in the noiseless case, the neurons responses can become synchronized. However, in the noisy case responses becomes more asynchronous, especially after the stimulus offset as the time approaches the expected reward time. The simulations shown in Figures 2A,B illustrate the ability of the network to represent different times. Longer decay times are implemented by stronger recurrent weights. However, for large enough weights the firing rates in the system no longer decay and it becomes bi-stable as it was shown in Gavornik and Shouval (2010). Although, for a deterministic system with infinitesimal precision it is possible to implement any decay time, for stochastic systems the proximity to the bi-stable regime imposes a practical upper limit to the reward time the system can represent (about 2.0–3.0 s).

The balancing of the LTP and LTD traces at the learned time is accomplished in part due to the difference in the characteristic time constants of the corresponding traces. The assumption that the LTP trace decays slower than the LTD trace i.e., τ^p > τ^d has experimental support (He et al., 2015). A second condition required is that the saturation of the LTD trace is higher than that of the LTP trace. This can be understood using the simple model illustrated in Figure 1A and the related equations. For instance, during the falling phase, the time dependence of the traces can be written as

where T~a (a ∈ {p, d}) is the saturation of the corresponding trace at t = t_stim (see Equation 4). If both traces balance at t = t_r then these two equations combine to give the relationship

From here it follows that since τ^p > τ^d it must be true that T~d>T~p, moreover as illustrated in Figure 2C, traces effectively saturate near their maximally allowed values Tmaxa, and thus it follows that Tmaxd>Tmaxp is necessary for convergence.

Recently, a different protocol has been used for training rodents, a protocol that allows the researcher to estimate when the animals expect to receive a reward. In this protocol, a visual stimulus is delivered, the animals decide when to act (lick) after the visual cue, and the longer they wait the more reward (water) they get, up to a maximal reward time t_max after which they receive no reward (Figure 3A, red line) (Namboodiri et al., 2015). Animals trained with this protocol learn to delay their licking such that it is close to t_max, and their chosen lick times are nearly optimal given the distribution of response times (Namboodiri et al., 2015). This protocol also results in changes to firing properties, similar to changes observed in previous protocols. Many of the cells in V1 exhibit a sustained increase in firing that terminates close to the time of action, and time of reward (Namboodiri et al., 2015).

We set out to examine if our two-trace learning rule is able to train a network with this training paradigm as well. To do this we used a similar network as in our previous simulations, and to this network we added a decision process that will trigger action. We used a very simple decision process such that once the average activity of excitatory cells in the network dips below a threshold (<15 Hz), action is initiated. We also assumed that reward is delayed by 200 ms, this delay is due to the time to initiate and execute the lick and also to the delay between the lick and the time that the reward signal arrives in the cortex (Hangya et al., 2015). While a delay is biologically realistic, some delay is also necessary for the simulations to succeed, and this delay must be larger than the crossing time between the LTP and LTD eligibility traces (see for example the crossing between the traces in Figure 2C, bottom).

With these assumptions, the network learns to increase synaptic efficacies and consequently increase the decay times of the neuronal responses. However, we find that this learning rule as it stands is insufficient to account for all the data. For instance, if the network is trained to one time interval such that the network initiates the action near the target time and then the reward time is changed to a shorter one, as illustrated by the upper and lower panels of Figure 3A, then the network will not be able to learn the new shorter reward time. This is because the action would be initiated at a time longer than the maximum reward time (lower panel in Figure 3A) and thus no reward would be given, and with the current learning rule the absence of reward results no change in synaptic efficacies. To address this we now consider, as exhibited experimentally (He et al., 2015), that there are two neuromodulators, one for expressing the LTP trace and the other for the LTD trace. We assume then that both of these neuromodulators are released when reward is delivered, however a much smaller release of the LTD associated neuromodulator is also assumed to be released upon action, even in the absence of reward.

Mathematically this means that if there is reward, we still use Equation (13) with Rp(t)=Rd(t)=R(t)δ(t-treward) at the time of reward. Where R(t) that defines the magnitude of reward has the form R(t) = R_max(t/t_max) for 0 < t < t_max and zero otherwise. However, if there is no reward, but action is initiated R^p(t) = 0 but Rd(t)=κRmaxδ(t-taction). The parameter κ << 1 scales the magnitude of neuromodulator in the absence of reward to be much smaller than that during reward, and this neuromodulator is delivered at the time the network initiates an action (t_action). Figure 3B illustrates the neuromodulator release profile. During the reward interval, both neuromodulators are delivered (cyan for LTP and magenta for LTD) but only R^d during the action (the scale used here is intended for clarity and differs in the actual simulations). Although, the assumption that a moderate amount of LTD related neuromodulator is released in the absence of reward was chosen for mathematical convenience, it is consistent with experimental observations. The role of an LTD related neuromodulator could be played by serotonergic neurons in the dorsal raphe nucleus (DRN) which innervate the visual cortex (Koh et al., 1991). These neurons have been shown to respond differentially to the delivery or omission of reward (Li et al., 2013, 2016). Moreover, as shown in Li et al. (2016) DRN neurons ramp-up during the waiting time of delivery of reward, regardless of whether the reward is delivered, however they respond in a phasic way when the reward is delivered, but slowly decay to baseline activity.

In Figure 4B we show the dynamics in a typical trial of the excitatory cells in this network after it is trained with the ramp-reward protocol with t_max = 1500 ms (red) (see Figure 4A). These dynamics can be altered when t_max is subsequently set to 1000 ms (Figure 4B, blue). Note that this reversal does not occur in a model in which only a single neuromodulator is used. The network dynamics change from trial to trial due to noise in the systems, consistent with experimental observations (Namboodiri et al., 2015). In Figure 4C we show summary statistics of the times that networks trained to 1500 ms (1) and 1000 ms (2) and cross the threshold from above. The medians obtained are 1308 ms for a network trained to 1500 and 875 ms for a network trained to 1000 ms. Note that we did not add the action delay to these numbers. These two distributions are statistically significantly different using an unpaired t-test, with P < 10⁻¹⁰.

The problem with the single neuromodulator model goes beyond its inability to transition back from long to short action times. In the absence of LTD expression, the response of a network trained from short to long action times will drift up. Eventually the network response time will move beyond t_max. We can see this by simulating a network with a t_max = 1000 ms, setting κ = 0.0 and running the simulations for many more iterations. In such a case (see Figure 4D, right) the network median threshold crossing drifts up until it reaches ~1090 ms. The action will be delayed with respect to the neural response, which means that the action will be observed at longer that 1200 ms, and beyond t_max. The slow drift to longer response times only ends when action times on every trial exceed t_max.

Experimentally the distribution of action times peaks close to the optimal time to obtain maximal reward, given the width of the distribution of action times (Namboodiri et al., 2015). One could ask if a network trained with CRF and two distinct neuromodulator accomplishes this as well. We find that properties of the action time distribution, in our model, depend on the parameters of the model, most specifically on the value of κ. In Figure 4D we show how the properties of the distribution depend on the value of κ. As κ decreases, the location of the peak of the distribution increases. Therefore, our model does not obtain optimal performance for all parameters, but could approach optimality for a specific value of κ.

The recurrent network described above can code for the expected time of reward, but we have not shown that it can code for the magnitude of reward. To investigate this, we set up a one layer network that is exposed to various 2D input patterns (of size N x N) (Figure 5A). Only one of these inputs is rewarded. The reward can be delivered either continuously or discretely through the time during which the pattern is presented. We examine two simple cases, one in which the reward is delivered at one point in time, and another in which it has a constant value throughout the presentation of the stimulus.

At first we examine the consequences of training with a discrete reward signal for expressing the LTP and LTD traces of the form: Rp(t)=Rd(t)=Rδ(t-treward). In this example the pattern is presented and remains on even after the time of the reward, which occurs at t_reward = 250 ms (see green line in Figures 5B,C). We find that this results in a potentiation of the weights that are stimulated by the rewarded pattern, while other weights are unchanged (Figure 5). Consequently the neuron becomes selective to the rewarded pattern. The strengthening of those synapses terminates once the LTP trace (T^p) and the LTD trace (T^d) are equal at the time of reward. Initially the weights vector is random (Figure 5B, left), the postsynaptic cell is only moderately activated (Figure 5B, center), and the LTP related trace T^p (solid line) is larger than the LTD trace T^d (dashed line) (Figure 5B, right), and both are much lower than saturation. These temporal profiles of the traces are determined by the model's parameters. For the T^p to increase faster than T^d we set η^p > η^d. This faster increase of T^p is essential for the initial growth of weights stimulated by the rewarded pattern. Learning causes the weight vector to resemble the rewarded pattern (Figure 5C, left), and the postsynaptic activation due to that pattern is significantly increased (Figure 5C, center). Note that when learning stabilizes the two eligibility traces at the time of reward have the same value (Figure 5C, right) and this results in no further changes (on average) to synaptic efficacies.

It can be useful for the response magnitude to learn to represent the expected reward magnitude. One could assume that the magnitude of the neuromodulators depends on the reward magnitude. If both R^p and R^d depend on reward magnitude in the same way, the fixed point will be independent on the reward magnitude because the steady state depends on the ratio between R^p and R^d. However, if they depend differently on reward magnitude then possibly reward magnitude could be represented in the response magnitude of the cells. We examine a simple case in which R^p depends on reward magnitude but the value of R^d is constant; it depends on the existence of reward, but not on the magnitude of reward.

To examine the impact of changing the value of R^p/R^d, we ran several simulations in which we modified this ratio. The simulations were run until the system reached a fixed point, and no further significant changes in the receptive fields were observed. As the value of R^p/R^d increased, so did the firing rate of the post-synaptic neuron in response to the rewarded pattern (see examples in Figure 6A), while the firing rates of non-rewarded patterns did not change significantly. On average the fixed point will occur when Rp(treward)Tijp(treward)=Rd(treward)Tijd(treward). Since the reward occurs during the rising phase of the traces, the fixed point (see Equation 8) is reached when:

Note that according to Equations (10, 11) T~p, T~d as well as τ~p and τ~d depend on H. Therefore, solving Equation (17) is equivalent as solving for H at the fixed point. The Hebbian term, in its simplest form, is proportional to the activity of the output neuron, therefore knowing the value of H gives an indirect measure of the network's response. Consequently, by solving numerically for H at different levels of R^p/R^d one can find how the magnitude of the Hebbian term (and thus the network response) at fixed point depends on the reward. A fixed point for this equation can be found if RpT~ijp<RdT~ijd while at the same time τ~p<τ~d. Although, T~p<Tmaxp and T~d<Tmaxd (see Equation 10), they typically saturate close to these maximal values, due to the large value of the Hebbian term with respect to the saturation levels, i.e., H≫Tmaxa. This implies new approximately necessary conditions for a fixed point RpTmaxp<RdTmaxd. Note that the value of τ~p and τ~d depends on τ^p and τ^d respectively, but also on η^p and η^d. By choosing appropriate values of η^p and η^d we can satisfy the condition τ~p<τ~d.

Although, Equation (17) cannot be solved in closed form, we can rewrite it in a simpler way by using the condition H≫Tmaxa to give

where r = R^p/R^d, k=Tmaxd/Tmaxp and ϕa=treward/(τaTmaxa). From this expression we can see that as r approaches k from below, then H → ∞, indicating that an increase in reward can be represented by the network as an increase in the network's response expressed in the Hebbian term. Solutions for r > k do not exist and this is indicated in Figure 6 by a sudden drop of the function to zero.

In Figure 6B we compare results of simulations to the solution of Equation (17), for Tmaxd/Tmaxp=1.5. We find a good agreement between the simulation results, averaged over many synapses, to the solution obtained by directly solving Equation (17) (solid line). We also find a large variability in each single trial between the value of H at different synapses (the error bars represent standard deviation over synapses), as well as between the value of H across trials within the same synapse. The variability primarily stems from the highly stochastic nature of presynaptic firing in each presynaptic neuron of each trial.

In Figure 6C we show the dependence of H on R^p/R^d for different values of Tmaxd/Tmaxp. These results confirm the stability condition (RdTmaxd>RpTmaxp) derived above.

In Figure 6A we showed via simulations how the postsynaptic firing rate depends on the R^p/R^d ratio, but our analysis in Equation (17) and in Figures 6B,C is for the Hebbian term H. Note that the fixed point for H is the same, independent of how H is actually calculated in simulations, as long as H is a form that is able to train the network and reach the fixed point. The details of the Hebbian term H do matter however for the post-synaptic rate. Given a specific functional form of H and a given presynaptic firing rate for the rewarded pattern, one could invert H to find the post-synaptic rate. For example, if H is simply H = c · r_pre · r_post, where r_pre and r_post are the pre- and post-synaptic firing rates, respectively and c a proportionality constant. Given that we define the fixed point as: H = H^fp the postsynaptic rate will take the form: rpost=Hfp/(c·rpre). Clearly for more complex Hebbian terms, for instance terms that depend on spike times, this inversion will be more complicated. However, the fixed point for H ≠ 0, will only be reached if the functional form of H chosen can reach this fixed point.

The results until now assumed a reward delivered at one point in time (a δ function reward). This makes analysis simpler but is clearly unrealistic. Another simple option is that reward is delivered uniformly throughout and extended period between times t₁ and t₂ such that t₂ − t₁ = Δt. These times represent the time interval during which reward is delivered, as illustrated in Figure 7A by the green bar at the bottom. Note that, as in the previous example, the reward is delivered while the pattern is still on. In that case the fixed point of Equation (13) takes the form:

What this means is that the fixed point is attained when the integral of the LTP trace times its reward magnitude is equal to the integral of the LTD trace times the LTD reward magnitude (see Figure 7B). From these assumptions, one derives a fixed point equation of the form:

This equation, like in the δ function reward case, is implicitly a function of H and can be solved for H. We do not have a closed form solution for this equation, but it can easily be solved numerically. Note that to obtain this expression we made the approximation that H is constant throughout the integral. This approximation is worse here than for the δ function case and therefore one should expect less agreement between the simulations and the solution of Equation (20). In Figure 7B we compare the solutions of this equation to simulation results. There is a surprisingly good agreement between simulations and analysis despite the rough approximations. Not surprisingly, the agreement is worse than for the δ function reward shown above. Importantly, and irrespective of the precise details, we find that it is possible to learn the magnitude of reward with two different neural modulators and this is possible for different temporal profiles of the reward signal.

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.