Brain-inspired classical conditioning model

Summary Classical conditioning plays a critical role in the learning process of biological brains, and many computational models have been built to reproduce the related classical experiments. However, these models can reproduce and explain only a limited range of typical phenomena in classical conditioning. Based on existing biological findings concerning classical conditioning, we build a brain-inspired classical conditioning (BICC) model. Compared with other computational models, our BICC model can reproduce as many as 15 classical experiments, explaining a broader set of findings than other models have, and offers better computational explainability for both the experimental phenomena and the biological mechanisms of classical conditioning. Finally, we validate our theoretical model on a humanoid robot in three classical conditioning experiments (acquisition, extinction, and reacquisition) and a speed generalization experiment, and the results show that our model is computationally feasible as a foundation for brain-inspired robot classical conditioning.


INTRODUCTION
Classical conditioning is regarded as a basic learning method for animals in which an association is built between a conditioned stimulus (CS) and a conditioned response (CR). The best-known experiment of classical conditioning was performed by Pavlov (1927). When a dog is presented with food (unconditioned stimulus, US), it will start to salivate (unconditioned response, UR). In Pavlov's research, a dog would hear a tone (CS) before being presented with food every time. After a number of trials, the dog started to salivate (CR) upon hearing a tone.

Computational model
Classical conditioning has attracted the interest of many researchers, and attempts have been made to build a computational model to reveal its mechanism. Rescorla and Wagner (1972) presented the first computational model of classical conditioning, which is named the Rescorla-Wanger model. This model can predict some important classical phenomena of classical conditioning, and it has led to a great deal of research, including modifications and alternative models. The Sutton-Barto (SB) model (Sutton and Barto, 1981) is a temporally refined extension of the Rescorla-Wanger model. This model learns to increase its response rate in anticipation of increased stimulation, and it can account for more phenomena observed in classical conditioning than the Rescorla-Wanger model can; furthermore, it has served as the precursor of many later models. The temporal difference (TD) model (Sutton and Barto, 1987) is an extension of the SB model. This model takes the form of a temporal difference prediction method and can successfully model the inter-stimulus interval (ISI) effect, which is regarded as a primary real-time effect of classical conditioning. Harry Klopf (1988) proposed a neuronal model that is modified from the Hebbian model to be more consistent with animal learning phenomena; this model can predict a wide range of classical conditioning phenomena. Schmajuk and DiCarlo (1992) presented a multilayer neural network called the S-D model. They mapped the nodes and connections onto regional cerebellar, cortical, and hippocampal circuits to obtain a model that can correctly describe the effects of hippocampal and cortical lesions on conditioning. Balkenius and Moren (1999) described a computational model of classical conditioning that is built on the assumption that the goal of learning is the prediction of a temporally discounted reward or punishment based on the current stimulus situation; notably, this model is well suited for robotic implementation. Johansson and Lansner (2002) presented an associative model of classical conditioning that is composed of a number of interconnected Bayesian confidence propagation neural networks (BCPNNs) implemented on the basis of Hebbian learning, and the output of this model fits the results of classical conditioning experiments. Zuo et al. (2005) introduced a spiking-neuron-based cognitive model. This model can simulate the learning process of classical conditioning with a reflex arc structure and a reinforcement learning method based on the Hebb rule, and an inverted pendulum experiment validated the effectiveness of this model.  built a model with classical conditioning behaviors based on a Bayesian network (CRMBBN). This model constructs cause-effect relationships between classical sensing and nonclassical conditional sensing by means of a Bayesian network, and it can successfully simulate many phenomena, such as acquisition, inter-stimulus effects, and extinction. Liu and Ding (2008) presented a dynamic policy adaptation framework inspired by classical conditioning. This model can successfully realize the self-learning process of classical conditioning and achieve adaptive network policy management. Antonietti et al. (2017) developed a detailed spiking cerebellar microcircuit model that can reproduce eyeblink classical conditioning and successfully fits real experimental datasets from humans.
Here, we build a brain-inspired classical conditioning (BICC) model that integrates and adopts existing biological findings of classical conditioning. With a quaternionic-rate-based synaptic learning rule, which is equivalent to spike-timing-dependent plasticity (STDP) ) on a timescale of seconds, the BICC model could predict the majority of classical conditioning phenomena. The computational model of biological classical conditioning enables a robot gets similar learning behavior and the capability of speed generalization. Furthermore, the changes in synaptic weights in this model may hint at the biological mechanism of classical conditioning.

Classical experiments
There are several classical conditioning experiments that can be used to verify the effectiveness of computational models. To enable reasonable comparisons with other well-known computational models, we follow the classical conditioning experiments outlined by Balkenius and Moren (1998). The stimulus before the + presents first, and the stimulus after the + presents later. The parenthesis indicates that the stimuli are presented and end simultaneously. The 0 indicates the result of training. The / indicates the result of the stimulus.

Acquisition
Acquisition is the ability to establish an association between a stimulus and a response, and it is the most fundamental process in classical conditioning. In an acquisition experiment, a CS is presented first and a US is presented subsequently after a small time interval for several trials; then, the response will be induced when the CS is presented on its own. This acquisition progress can be described as follows. For an acquisition experiment involving an eyelid response in the albino rabbit, the response level forms an S-shaped curve similar to a sigmoid function (Balkenius and Moren, 1998;Schneiderman et al., 1962).

Inter-stimulus interval effect
The ISI effect is a primary real-time effect of classical conditioning (Sutton, 1990). The ISI represents the time interval between the presentation of the CS and US, and it can be divided into three types (Balkenius and Moren 1999): delay conditioning A, delay conditioning B, and trace conditioning ( Figure S1). In delay conditioning A, the US is presented immediately when the CS terminates. In delay conditioning B, the CS is still present when the US is presented, and the CS and US terminate simultaneously. In trace conditioning, the CS and US have fixed lengths, and the CS terminates before the onset of the US. In empirical studies conducted by Schneiderman and Gormezano (1964) and Smith et al. (1969), the ISI-CR frequency function revealed a concave-down shape during the acquisition and extinction process.

Extinction
In an extinction experiment, the acquired response will disappear gradually if only a CS is presented without the subsequent US. The extinction process can be described as follows:

Reacquisition effect
The reacquisition effect is demonstrated when an animal relearns a previously extinguished association, and the relearning phase is faster than the initial learning phase. iScience Article Blocking Blocking refers to the following phenomenon: when a first stimulus CS 1 has been associated with a response and a second stimulus CS 2 then is presented and ends simultaneously with CS 1 , the attempt to associate the second stimulus CS 2 with the response will be blocked. Blocking experiments show that the association of a stimulus with a response is not independent of earlier learning. The blocking process can be described as follows, where the parentheses are used to indicate that CS 1 and CS 2 are presented and end simultaneously.

Secondary conditioning
In a secondary conditioning experiment, CS 1 has been associated with the response induced by the US, and CS 1 is then used as the US for CS 2 to build an association to the response. The effect of such secondary conditioning is typically weak, and CS 1 will undergo extinction, whereas CS 2 is reinforced. The secondary conditioning progress can be described as follows:

Conditioned inhibition
In a conditioned inhibition experiment, CS 2 and CS 0 have been associated with the response, and then a third stimulus CS 1 is presented and ends simultaneously with one of the previously conditioned stimuli CS 0 without the US. In the test phase, CS 1 will inhibit the ability of CS 2 to induce the response. The conditioned inhibition process can be described as follows, where the parentheses are used to indicate that the stimuli are presented and end simultaneously.

Facilitation by an intermittent stimulus
Under normal acquisition conditions, the conditioning to CS 1 is weak in the case of trace conditioning due to the long ISI. Under conditions of facilitation, the conditioning to CS 1 is facilitated by an additional CS 2 . The facilitation process can be described as follows: Normal The CS 1 and CS 2 are presented and ended simultaneously under the conditions of overshadowing; the associative strength acquired by CS 1 and CS 2 are weaker than the CS 1 or CS 2 conditioned alone in the normal acquisition condition (Angulo et al., 2020). The overshadowing process can be described as follows: Normal The CS 1 and CS 2 have been associated with the response, respectively, then the following CS 1 -CS2 presentations result in a weight decrement (Rescorla and Wagner, 1972). The overexpectation process can be described as follows:

Recovery from overshadowing
In the overshadowing experiment, the extinction of the CS 1 will lead to an increased responding to the CS 2 (Matzel et al., 1985). The recovery from overshadowing process can be described as follows:

Recovery from forward blocking
In the forward blocking experiment, the extinction of the blocker CS 1 will lead to an increased response to the blocked CS 2 (PinenO et al., 2005). The recovery from forward blocking process can be described as follows:

RESULTS
The neural circuit underlying delay eyeblink conditioning has been well described in (Hansel et al., 2001;Ten Brinke et al., 2019;Wang et al., 2018;Hogri et al., 2015;Takehara-Nishiuchi, 2018), and we propose our BICC model based on these findings. The architecture of the BICC model is shown in Figure 1.

Model evaluation
In this section, we use the changes in synaptic weight between the PN (pontine nuclei) and the IPN (interpositus nucleus) to evaluate this model; NP CSi represents the neuron population of the corresponding CSin the PN, and NP R represents the neuron population that controls the response in the IPN. The PN deliver the information from the CS, and the IPN generates the CR via the motor control pathway, as introduced in the Methods section.

Inter-stimulus interval effects
We use both delay and trace conditioning experiments to test our model, and the variation in the synaptic weight between NP CS and NP R is shown in Figure 2A. The curves initially show a marked increase and then exhibit a concave-down shape, which is consistent with the results of the rabbit experiment presented in Schneiderman and Gormezano (1964) and Smith et al. (1969). There is an optimal interval time for learning under every condition in this model, being 1.3 s for trace conditioning and 3.1 s for both types of delay conditioning. iScience Article Learning curves Figure 2B shows the results for the learning curves in the model. Under each condition, the curve is an S-shaped acquisition curve. With fewer than eight trials, the synaptic weight between NP CS and NP R is greater in the case of delay conditioning B than in the case of delay conditioning A, and the smallest weight is observed in the trace conditioning case. As the number of trials increases, the synaptic weight increases to a stable constant; ultimately, the synaptic weight under delay conditioning A is greater than that under trace conditioning, and the smallest final weight is observed in the case of delay conditioning B. We therefore select delay conditioning A for the remainder of the experiments unless otherwise stated, with both the CS and the US continuing for 2 s. Figure 2C shows the variations in synaptic weight between NP CS and NP R in the acquisition and extinction experiments. In the acquisition experiment, the CS is presented at 0 s and ends at 2 s and the US is presented at 2 s and ends at 4 s. In the extinction experiment, only a CS is presented at 0 s and ends at 2 s, without a US. The results of the reacquisition experiment are shown in Figure 2D. It is easy to see that in the reacquisition experiment, fewer (1) In the delay conditioning A, the duration of the CS varies from 0 to 6 s, and then the US is presented immediately after the CS and continues for 2 s. The ISI is equal to the length of the CS.

Acquisition, extinction, and reacquisition experiments
(2) In the delay conditioning B, the duration of the CS varies from 1 to 7 s, whereas the US continues for 1 s, and the CS and US end simultaneously. The ISI is equal to the difference between the length of the CS and the length of the US.
(3) In trace conditioning, the CS and US each continue for 2 s. The ISI is equal to the difference between the start time of the CS and the start time of the US, and it varies from 0 to 6 s. (B) The learning curves in the model. (1) In delay conditioning A, the CS and US each continue for 2 s, so the interval time is 2 s. (2) In delay conditioning B, the CS continues for 3 s and the US continues for 1 s, so the interval time is 2 s. (3) In trace conditioning, the CS is presented at 0 s and continues for 2 s, and the US is presented at 2.5 s and continues for 2 s, so the interval time is 2.5 s.
(C) The variations in synaptic weight between NP CS and NP R . W Acquisition1 and W Acquisition2 represent the weight variations in the acquisition experiment when the initial weight is 0 and 0.5, respectively. W Extinction represents the weight variations in the extinction experiment, where there is only a CS.
(D) Reacquisition experiment. It is easy to see that fewer trials are needed in the reacquisition condition (fewer than 10 trials) than in the acquisition condition (approximately 25 trials) to achieve the same weight.

OPEN ACCESS
iScience 24, 101980, January 22, 2021 5 iScience Article trials are required to achieve a higher synaptic weight between NP CS and NP R than in the acquisition experiment. This is because during acquisition or extinction, not only the synaptic weight but also the number of synapses involved changes. In the reacquisition experiment, more synapses are involved in the synaptic weight changes, so the learning rate is faster than in the acquisition stage. We combine the acquisition, extinction, and reacquisition experiments into a single overall experiment. The results are shown in Figure S2. Figure 3A shows the results of the blocking experiment. It is easy to see that the synaptic weight between NP CS2 and NP R is too small to induce a response. In the blocking stage, the changes in single synaptic weights between NP CS1 and NP R and between NP CS2 and NP R are identical because of the synchronization of CS 1 and CS 2 , but there are many more synapses involved in NP CS1 than in NP CS2 , which causes the change in W R;CS1 to be much greater than that in W R;CS2 .

Blocking experiment
Secondary conditioning experiment Figure 3B shows the results of the secondary conditioning experiment. Here, CS 1 is treated as the US to build an association between CS 2 and the response, and the corresponding synaptic weight is typically weak because the synaptic weight between NP CS1 and NP R exhibits an extinction effect at the same time. In the first 10 trials, only CS 1 and US are presented to build a conditioned response. In the remainder of the trials, CS 1 and CS 2 are presented and end simultaneously, and then the US is presented subsequently. (B) Secondary conditioning experiment. In the first 25 trials, CS 1 is presented first and US is presented subsequently, as in the normal acquisition experiment. In the last 25 trials, CS 1 is treated as the unconditioned stimulus and is presented after CS 2 .
(C) Conditioned inhibition experiment. In the first 25 and second set of 25 trials, CS 2 and CS 0 , respectively, are combined with US for conditioning on the response. In the subsequent 25 trials, CS 0 and CS 1 are presented and end simultaneously, without US. In the last 25 trials, only CS 1 and CS 2 are presented. (D) Facilitation experiment. CS 1 continues for 2 s, CS 2 continues for 3 s, and the US continues for 1 s. During normal acquisition, CS 1 is presented first, and 2 s later, the US is presented. During facilitated acquisition, CS 1 is presented first, CS 2 is presented immediately when CS 1 ends, and 2 s after CS 1 ends, the US is presented. CS 2 and the US end simultaneously. Conditioned inhibition experiment Figure 3C shows the results of the conditioned inhibition experiment. In the first and second sets of 25 trials, excitatory synapse connections to NP R are built for NP CS2 and NP CS0 , respectively. In the third set of 25 trials, NP CS0 exhibits an extinction effect, whereasNP CS1 builds inhibitory synapse connections to NP R because of synchronism. In the last 25 trials, the number of synapses between NP CS1 and NP R increase because of the inhibitory connections and the negative weight changes. At the beginning of the last set of 25 trials, the inhibition effect from CS 1 is not sufficient to completely inhibit the excitatory input from CS 2 , so CS 2 can still induce the response. With the enhancement of the inhibition effect from CS 1 , at the end of the last set of 25 trials, CS 1 can inhibit the response induced by CS 2 . In the extinction experiment, more than 50 trials are needed for CS 2 to lose the ability to induce the response, whereas in the conditioned inhibition experiment, fewer than 25 trials are needed because of the inhibition effect from CS 1 . Figure 3D shows the results of the facilitation experiment. It is easy to see that the synaptic weight between NP CS1 and NP R is stronger under facilitated acquisition than under normal acquisition. Under normal acquisition, the synaptic weight is weak because of the long ISI. Under facilitated acquisition, CS 2 is conditioned on the response, and the response will be induced twice, by CS 2 and the US, thus causing CS 1 to build stronger synaptic connections.

Facilitation experiment
Overshadowing experiment Figure 4A shows the results of the overshadowing experiment. The synaptic weights between NP CS1 and NP R and NP CS2 and NP R are identical because of the synchronization of CS 1 and CS 2 . It is easy to see that the synaptic weight in the overshadowing condition is weaker than that in the normal acquisition condition. In the overshadowing condition, the NP R is stimulated by CS 1 and CS 2 simultaneously, and they contribute equally to build an association to response. So, the synaptic weight in the overshadowing condition is weaker, and it is about half of that in the normal acquisition condition.
Overexpectation experiment Figure 4B shows the results of the overexpectation experiment. The synaptic weight between NP CS1 and NP R and NP CS2 and NP R have decreased in the last 25 trials, and it provided that the following CS 1 -CS2 presentations result in a weight decrement. In the first and subsequent 25 trials, the CS 1 and CS 2 build a CR with US, respectively. In the last 25 trials, the CS 1 and CS 2 stimulate the NP R simultaneously. So, the firing rate of the response neuron increases faster and lasts longer, and it makes an extinction effect until the model is stable again.
Recovery from overshadowing experiment Figure 4C shows the results of the recovery from overshadowing experiment. The synaptic weight between NP CS2 and NP R has increased when the CS 1 ended, and it provided that the extinction of the CS 1 will lead to an increased respondse to the CS 2 . The NP R is only stimulated by CS 2 when the CS 1 ended, and there is an acquisition effect when the US is presented.
Recovery from forward blocking experiment Figure 4D shows the results of the recovery from forward blocking experiment. The synaptic weight between NP CS2 and NP R has increased when the CS 1 ended, and it provided that the extinction of the blocker CS 1 will lead to an increased response to the blocked CS 2 . Similar to the recovery from overshadowing experiment, the NP R is only stimulated by CS 2 when the CS 1 ended, and there is an acquisition effect when the US is presented.
The results of the presented model in the various experiments and the comparison results with existing models are summarized in Table 1.

Robotic classical conditioning experiments
We use acquisition, extinction, and reacquisition experiments and speed generalization experiment to evaluate this model on a humanoid robot.

Acquisition, extinction, and reacquisition experiments
We selected three classical conditioning experiments-acquisition, extinction, and reacquisition-to evaluate this model on a humanoid robot. A red robot was used as the participant robot, and the various stimuli used in the experiments are shown in Figure 5A. We use template matching to identify different stimuli. The sequences of stimuli in these experiments are shown in Figure 5B. In the acquisition experiment, the participant robot first was shown a red fist toy (CS), then was shown a blue robot (US), and subsequently took an avoidance action (UR). After learning, the participant robot would take an avoidance action (CR) when it saw the red fist toy (CS).
In the extinction experiment, only the CSwas presented; after several trials, the participant robot did not perform the CR upon seeing the red fist toy (CS). In the reacquisition experiment, the participant robot could establish the CR in fewer trials than in the acquisition experiment. The experimental results are shown in Figures 5C and 5D.

Speed generalization experiment
The speed generalization experiment is that the robot trained at a slow speed on a navigation task; then it could navigate the track at a higher speed without training (Herreros et al., 2013). The humanoid robot Nao is a biped (B) Overexpectation experiment. In the first 25 trials, the CS 1 is presented first and the US is presented subsequently. In the subsequent 25 trials, the CS 2 is presented first and the US is presented subsequently. In the last 25 trials, the CS 1 and CS 2 are presented and end simultaneously, and then the US is presented subsequently. (C) Recovery from overshadowing experiment. The first 25 trials are the overshadowing process, the CS 1 and CS 2 are presented and end simultaneously, and then the US is presented subsequently. The last 25 trials are the recovery process, the CS 1 is ended, and the CS 2 presented first and the US is presented subsequently. (D) Recovery from forward blocking experiment. The first 30 trials are the blocking process. In the first 10 trials, the CS 1 is presented first and the US is presented subsequently. In the subsequent 20 trials, the CS 1 and CS 2 are presented and end simultaneously, and then the US is presented subsequently. The last 20 trials are the recovery process, wherein the CS 1 is ended, and the CS 2 presented first and the US is presented subsequently. iScience Article robot. According to its documentation, the speed parameter is the fraction of the maximum speed, such as 100% means the full speed. According to our verification, may be due to the aging problem, the robot cannot move accurately at a given speed, and often deviates from the direction when walking in a straight line. So, we test the speed generalization capabilities of the model in both real and simulation environments.
The real environment is shown as Figure 6A. In the real environment, the robot is trained at 50% speed and tested at 100% speed. We use simple image recognition algorithms to recognize the CS (black line, Figure 6B) and the US (red line, Figure 6C). For the black line recognition, we retain the central area of the image, then convert the image to a binary image according to the given threshold value, delete the small-area object, and finally perform horizontal line detection to complete the recognition. For the red line recognition, we extract the matching color area by setting the threshold value and then detect the horizontal line to complete the recognition. The result of the real environment is shown as Figure 7.
The simulation environment is shown as Figure 8A. In the simulation environment, we use two experiments to test the speed generalization capacity of the model. In the first experiment, we test the induced time of CR and UR at a given speed in the range of 100%-200% with 10% speed increments. The result is shown as Figure 8B. In the second experiment, we test the induced distance of CR at a given speed at 100%, 120%, 150%, 200%, 300%, 350%, and 400%. The result is shown as Figure 8C. And if the speed is greater than 350%, the robot experiment fails.
The results of real and simulation environments show that the proposed model has the capacity of speed generalization.

DISCUSSION
The BICC model exhibits long-term depression (LTD) at GC (granule cell)-PU (Purkinje cell) synapses and long-term potentiation (LTP) at PN-IPN synapses, which is consistent with the findings from electrophysiological experiments on eyeblink conditioning presented in Koekkoek et al. (2003) and Steuber et al.
This table is adapted from Johansson and Lansner (2002) and Balkenius and Moren (1998). An * means that the model could reproduce the correlation feature, o means that the model could reproduce partially, and -means that it is not mentioned or was unable to reproduce. SB (Sutton and Barto 1981)

OPEN ACCESS
iScience 24, 101980, January 22, 2021 9 iScience Article (2007) and Raman, (2006, 2008), respectively. LTD is also exhibited at PN-IPN synapses in this model, representing another mechanism of synaptic plasticity that is consistent with the electrophysiological experiments reported in Zhang and Linden (2006), but there have been fewer reports of this mechanism than of the former.

LTD at GC-PU synapses
If only the CS is presented before learning, the GC (granule cell) receives the CS from the PN and then projects it to an Int.N (inhibitory interneuron). The inhibitory input from the Int.N will inhibit the excitatory input from the GC and the spontaneous firing of the PU, so the firing rate of the PU will immediately drop to zero; in other words, the firing of the PU will be paused. The synaptic weight between the GC and PU will not change because the postsynaptic neurons in the PU is not fired.
If only the US is presented before learning, the US is projected to the motor control pathway and induces the UR. In parallel, the IO (inferior olive) receives the US from the US pathway and then projects it to the IPN and PU via the cf (climbing fiber). The excitatory input from the IO to the PU will enhance the inhibition from the PU to the IPN. The synaptic weight between the GC and PU remains unchanged because the presynaptic neuron in the GC is not fired.
In the acquisition experiment, the firing rate of the PU will gradually decay to zero because of the additional excitatory input from the US. With the increased firing rate of the presynaptic neurons in the GC and the iScience Article decreasing firing rate of the postsynaptic neurons in the PU, the synaptic weight between the GC and PU will decrease and exhibit LTD. An intuitive explanation of how spike-based STDP can influence synaptic efficiency through a rate-based mechanism can be found in Bengio et al. (2017). In short, the synaptic weight is updated in proportion to the product of the presynaptic firing rate and the temporal rate of change in activity on the postsynaptic side. The synaptic weight is updated in our model based on more factors, as well, as detailed in . The synaptic weight between the GC and PU will decrease with repeated pairings of the CS and US, whereas the inhibitory input from the Int.N and the excitatory input from the IO remain unchanged, so the Purkinje cell will pause spontaneous firing. This phenomenon is consistent with the electrophysiological experiments on eyeblink conditioning presented in Wetmore et al. (2014) and Hansel et al. (2001).

LTP and LTD at PN-IPN synapses
The single synaptic weight changes in a single trial in the acquisition and extinction experiments are shown in Figure 9A, and the firing rates of the neurons in the acquisition and extinction experiments are shown in Figures 9B and 9C, respectively. In the acquisition experiment, the change in the synaptic weight is negative because the temporal rate of change of the postsynaptic neurons in the IPN, which is represented as V R in Figures 9B, is smaller than that of the presynaptic neuron in the PN, which is represented as V CS in Figure 9B, from 0 to 2 s. The CS ends at 2 s; then the US is presented, continues for 2 s, and finally ends at 4 s. From 2 to 4 s, the firing rate of the presynaptic neuron decreases because the CS has ended, whereas the firing rate of the postsynaptic neuron increases because the US is being presented. The change in the synaptic weight is positive because of the decreasing firing rate of the presynaptic neuron and the increased firing rate of the postsynaptic neuron. This model exhibits the acquisition effect if the positive term is greater than the negative term and exhibits the extinction effect if the positive term is less than the negative term, and the model reaches a steady state when the positive term is equal to the negative term.
In the acquisition experiment, the number of excitatory synapses between the PN and IPN increases, which is consistent with the electrophysiological experiments on eyeblink conditioning presented in Kleim et al. (2002) and Weeks et al. (2007).
Our model suggests that the cerebellar cortex, especially the IPN, plays a critical role in classical conditioning. In our model, the LTD at GC-PU synapses leads to a reduction in the excitatory input from the GC. Although there is an excitatory input from the IO when the US is presented, the PU will be paused due to the loss of excitatory input from the GC. In the BICC model, classical conditioning can be achieved without the PU but not without the IPN, which is consistent with Lavond and Steinmetz (1989).

Limitations of the study
Our model cannot reproduce the experiment of latent inhibition (Lubow and Moore, 1959), spontaneous recovery (Pavlov, 1927), and unblocking (Bradfield and Mcnally, 2008). The latent inhibition effect is that a familiar stimulus takes longer to build an association to response than a new stimulus. Our model cannot reproduce the latent inhibition experiment because the model doesn't distinguish between familiar stimulus and new stimulus. The spontaneous recovery effect is the reappearance of a response that had been extinguished. Our model cannot reproduce this experiment, and we think that the spontaneous recovery effect may require involvement of more brain regions. The unblocking effect is that the responding of the blocked CS 2 increases by increase in the intensity or duration of the US, or increase in the number of the US. Our model failed in the unblocing effect experiment, if the CS 1 has built an association to response, the CS 2 cannot build the association, no matter how the US changes. In the future, we will improve our model to reproduce more experiments.

Resource availability Lead contact
Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Yi Zeng (yi.zeng@ia.ac.cn).

Computational Model
The architecture of the BICC model is shown in Figure 1. The change in the firing rate of a neuron that has received a conditioned stimulus (CS) or an unconditioned stimulus (U S) in this model can be described as shown in Equation 1. S j (t) represents stimulus j ∈ {CS 1 , CS 2 , ..., CS N , US} at time t, and V j (t) represents the firing rate of the neuron. We set S j (t) to 1 if the stimulus is present at time t and set it to 0 if the stimulus has ended. The firing rate of the neuron will increase when the stimulus is present and decay when the stimulus has ended. C is the parameter that controls the rate of increase and decrease.
Synaptic plasticity occurs in two areas: between GC (granule cell) and PU (Purkinje cell) and between the PN (pontine nuclei) and IPN (interpositus nucleus). The Purkinje cell is the sole output of the cerebellar cortex and inhibits the neurons in the IPN (interpositus nucleus) (Takehara-Nishiuchi, 2018), and it shows spontaneous firing due to its intrinsic membrane properties, as supported by evidence from several studies (Häusser and Clark, 1997;Raman and Bean, 1997). In the cerebellar cortex, the firing rate of the Purkinje cell depends on its own spontaneous firing, the stimulation from the CS and U S, and the inhibitory information from inhibitory interneurons, as described in Equation 2.
is the firing rate of the Purkinje cell at time t, and V pu−self is the spontaneous firing rate of the Purkinje cell itself. For k ∈ {CS 1 , CS 2 , ..., CS N }, W gc k (t) is the synaptic weight between the GC (granule cell) and the PU (Purkinje cell), V gc k is the firing rate of the GC (granule cell), W int.n is the synaptic weight between the Int.N (inhibitory interneuron) and the PU (Purkinje cell), and V int.n k is the firing rate of the Int.N (inhibitory interneuron). W io is the synaptic weight between the IO (inferior olive) and the PU (Purkinje cell), and V io (t) is the firing rate of the IO (inferior olive). W int.n and W io remain unchanged in this model.
In the IPN (interpositus nucleus), the firing rate of the neuron depends on the inhibitory information from the Purkinje cell and the stimulation from the CS and U S, as described in Equation 3.
is the firing rate of the neuron in the IPN (interpositus nucleus) at time t. W pu is the synaptic weight between the IPN (interpositus nucleus) and the PU (Purkinje cell), and V pu (t) is the firing rate of the Purkinje cell. W pn k (t) is the synaptic weight between the PN (pontine nuclei) and the IPN (interpositus nucleus), and V pn k (t) is the firing rate of the PN (pontine nuclei). W io is the synaptic weight between the IO (inferior olive) and the IPN (interpositus nucleus), and V io (t) is the firing rate of IO (inferior olive). W pu and W io remain unchanged in this model.
The synaptic plasticity in this model is defined as shown in Equation 4. W ij in this model represents the synaptic weight between postsynaptic neuron i and presynaptic neuron population j. The postsynaptic neurons contain the neurons in the PU (Purkinje cell) and the IPN (interpositus nucleus), and the presynaptic neurons contain the neurons in the GC (granule cell) and the PN (pontine nuclei). where K ij (T ) in this model represents the number of synapses between postsynaptic neuron i and presynaptic neuron population j. ∆w ij (t) is the incremental weight of a single synapse, which is calculated using Equation 5.
is the STDP function reported in Song et al., 2000;Shouval et al., 2002). And the general shape of the STDP function is shown as Figure S3.
Equation 5 is a quaternionic-rate-based synaptic learning rule derived from spike-timing-dependent plasticity, which we proposed in . It is easy to see that there are three main terms that affect the changes in synaptic weight, namely, the product of the pre-and postsynaptic neuron firing rates, the product of the presynaptic neuron firing rate and the time derivative of the postsynaptic neuron firing rate, and the product of the postsynaptic neuron firing rate and the time derivative of the presynaptic neuron firing rate, with corresponding coefficients α, β and γ. Equation 5 is used to eliminate the incompatibility of spike-timing-dependent plasticity and classical conditioning on a temporal scale, because the effective time window in spike-timing-dependent plasticity is on the millisecond scale, such as 40 ms, but it can be as long as several seconds in classical conditioning.
In the learning process of a biological system, changes occur not only in the form of synaptic weight changes but also in the forms of synaptic growth and synaptic elimination (Johnston et al., 2009;Bailey et al., 2004). The change in the number of synapses, ∆K ij (T + 1), is defined in Equation 6.
If W ij (T ) > 0, this means that most of the synapses between neuron population j and neuron population i are excitatory connections. If ∆w ij (T ) > 0, this means that the strength of excitatory connections increases and the number of excitatory synapses simultaneously increases, resulting in synaptic growth. If ∆w ij (T ) < 0, this means that the strength of excitatory connections decreases and the number of excitatory synapses simultaneously decreases, resulting in synaptic elimination.
If W ij (T ) < 0, this means that most of the synapses between neuron population j and neuron population i are inhibitory connections. If ∆w ij (T ) > 0, this means that the strength of inhibitory connections decreases and the number of inhibitory synapses simultaneously decreases, resulting in synaptic elimination. If ∆w ij (T ) < 0, this means that the strength of inhibitory connections increases and the number of inhibitory synapses simultaneously increases, resulting in synaptic growth.

Parameter Setting
We make the stimulus continues for 2s, and the characteristic behaviors of a neuron for which we set the parameter C to 1/200, 2/200, 3/200, 4/200 in this model is shown in Figure S4A. When the stimulus ends, the slow decrease of the neuron's firing rate is conducive to the establish the conditioned response in trace condition, but if it takes too long for the neuron's firing rate to decay to zero, it is not conducive to the convergence of the model, such as the parameter is 1/200. The neuron's characteristic behaviors are similar when the parameter is 2/200, 3/200, 4/200, and each parameter could make the model establish the association between the stimulus and the response. So we set the parameter C to 3/200 in this model which makes the firing rate of the neuron increase rapidly and decrease slowly.
The parameters α, β and γ calculated by the STDP function that reported in  can't make the model reproduce the classical experiments. So, we calculate these parameters using the incremental weights in three states.
We convert the Equation 1 to a function with t as the variable, and the converted function is shown in Equation 7. The t c is the continues time of the stimulus.
We simplified the equation for calculating the firing rate of postsynaptic neuron i, as shown in Equation 8.
We calculate the parameters α, β and γ using the incremental weights in three states, they are the acquisition state, the acquisition stable state and the extinction state. In the acquisition state and the acquisition stable state, the CS is presented at 0s and ends at 2s, and the U S is presented at 2s and ends at 4s. In the extinction state, only a CS is presented at 0s and ends at 2s, without a U S.
Because of the non-linear transformation in Equation 8, we calculate the firing rate of the postsynaptic neuron i in three states respectively, and then use the fourth-order Fourier function to fit the firing rate curve. So, we can get a function of postsynaptic neuron's firing rate with t as the variable, as shown in Equation 9. And the parameters in three states are shown in Table S1.
We set the ∆W ij (T ) to 0.25, 0, -0.1 in the acquisition state, the acquisition stable state and the extinction state respectively. By solving the equations, we can get the α = −0.0035, β = 0.3441 and γ = −0.5510. Consider that the fourth-order Fourier function cannot fully fit the firing rate curve of the postsynaptic neuron, we select parameters of α = −0.0035, β = 0.35, γ = −0.55 in this model. And the STDP function calculated according to these parameters is shown in Figure S4b.
When the W ij (T ) and ∆w ij (T ) are constant, the k determines the rate of synaptic growth and elimination. The larger the k, the faster the synaptic growth and extinction in the process of acquisition and extinction, and the model needs less trials to complete the acquisition or extinction. The smaller k is, the slower the rate of synaptic growth and extinction in the process of acquisition and extinction, and more trials are required for the model to complete acquisition or extinction. We select parameters of k = 8, so that the model can complete acquisition or regression with a relatively appropriate number of trials (such as 25), and facilitate observation of experimental phenomena.
Shouval, H. Z., Bear, M. F. and Cooper, L. N. (2002), 'A unified model of nmda receptor-dependent bidirectional synaptic plasticity', Proceedings of the National Academy of Sciences of the United States of America 99 (16)