Motor system model with reinforcement learning drives virtual arm (Dura-Bernal et al 2017)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:194897
"We implemented a model of the motor system with the following components: dorsal premotor cortex (PMd), primary motor cortex (M1), spinal cord and musculoskeletal arm (Figure 1). PMd modulated M1 to select the target to reach, M1 excited the descending spinal cord neurons that drove the arm muscles, and received arm proprioceptive feedback (information about the arm position) via the ascending spinal cord neurons. The large-scale model of M1 consisted of 6,208 spiking Izhikevich model neurons [37] of four types: regular-firing and bursting pyramidal neurons, and fast-spiking and low-threshold-spiking interneurons. These were distributed across cortical layers 2/3, 5A, 5B and 6, with cell properties, proportions, locations, connectivity, weights and delays drawn primarily from mammalian experimental data [38], [39], and described in detail in previous work [29]. The network included 486,491 connections, with synapses modeling properties of four different receptors ..."
Reference:
1 . Dura-Bernal S, Neymotin SA, Kerr CC, Sivagnanam S, Majumdar A, Francis JT, Lytton WW (2017) Evolutionary algorithm optimization of biological learning parameters in a biomimetic neuroprosthesis. IBM Journal of Research and Development (Computational Neuroscience special issue) 61(2/3):6:1-6:14
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Abstract Izhikevich neuron;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; GabaB; NMDA; AMPA;
Gene(s):
Transmitter(s): Glutamate; Gaba;
Simulation Environment: NEURON; Python;
Model Concept(s): Learning; Reinforcement Learning; Reward-modulated STDP; STDP; Motor control; Sensory processing;
Implementer(s): Dura-Bernal, Salvador [salvadordura at gmail.com]; Kerr, Cliff [cliffk at neurosim.downstate.edu];
Search NeuronDB for information about:  GabaA; GabaB; AMPA; NMDA; Gaba; Glutamate;
"""
MSARM

Code to connect a virtual arm to the M1 model
Multipe arm options depending on self.type:
- 'randomInput': the object provides random values matching the format (used for test/debugging)
- 'dummyArm': simple kinematic arm implemented in python 
- 'musculoskeletal': realistic musculoskeletal arm implemented in C++ and interfaced via pipes

Adapted from arm.hoc in arm2dms


Version: 2015jan28 by salvadordura@gmail.com
"""

from neuron import h
import arminterface
from numpy import array, zeros, pi, ones, cos, sin, mean
from pylab import concatenate, figure, show, ion, ioff, pause,xlabel, ylabel, plot, Circle, sqrt, arctan, arctan2, close
from copy import copy
from random import uniform, seed, sample, randint


[SH,EL] = [X,Y] = [0,1]
[SH_EXT, SH_FLEX, EL_EXT, EL_FLEX] = [0,1,2,3]

  
class Arm:
    #%% init
    def __init__(self, type, anim, graphs): # initialize variables
        self.type = type # 'randomOutput', 'dummyArm', 'musculoskeletal'
        self.anim = anim # whether to show arm animation or not
        self.graphs = graphs # whether to show graphs at the end


    ################################
    ### SUPPORT METHODS
    ################################

    # convert cartesian position to joint angles
    def pos2angles(self, armPos, armLen):
        x = armPos[0]
        y = armPos[1]
        l1 = armLen[0]
        l2 = armLen[1]
        elang = abs(2*arctan(sqrt(((l1 + l2)**2 - (x**2 + y**2))/((x**2 + y**2) - (l1 - l2)**2)))); 
        phi = arctan2(y,x); 
        psi = arctan2(l2 * sin(elang), l1 + (l2 * cos(elang)));
        shang = phi - psi; 
        return [shang,elang]

    # convert joint angles to cartesian position
    def angles2pos(self, armAng, armLen):        
        elbowPosx = armLen[0] * cos(armAng[0]) # end of elbow
        elbowPosy = armLen[0] * sin(armAng[0])
        wristPosx = elbowPosx + armLen[1] * cos(+armAng[0]+armAng[1]) # wrist=arm position
        wristPosy = elbowPosy + armLen[1] * sin(+armAng[0]+armAng[1])
        return [wristPosx,wristPosy] 

    #%% setTargetByID
    def setTargetByID(self, id, startAng, targetDist, armLen):
        startPos = self.angles2pos(startAng, armLen)
        if id == 0:
            targetPos = [startPos[0]+0.15, startPos[1]+0]
        elif id == 1:
            targetPos = [startPos[0]-0.15, startPos[1]+0]
        elif id == 2:
            targetPos = [startPos[0]+0, startPos[1]+0.15]
        elif id == 3:
            targetPos = [startPos[0]+0, startPos[1]-0.15]
        return targetPos

    #%% setupDummyArm
    def setupDummyArm(self):
        if self.anim:
            ion()
            self.fig = figure() # create figure
            l = 1.1*sum(self.armLen)
            self.ax = self.fig.add_subplot(111, autoscale_on=False, xlim=(-l/2, +l), ylim=(-l/2, +l)) # create subplot
            self.ax.grid()
            self.line, = self.ax.plot([], [], 'o-', lw=2)
            self.circle = Circle((0,0),0.04, color='g', fill=False)
            self.ax.add_artist(self.circle)

    #%% runDummyArm: update position and velocity based on motor commands; and plot
    def runDummyArm(self, dataReceived):
        friction = 0.5 # friction coefficient
        shang = (self.ang[SH] + self.angVel[SH] * self.interval/1000) #% update shoulder angle
        elang = (self.ang[EL] + self.angVel[EL] * self.interval/1000) #% update elbow angle
        if shang<self.minPval: shang = self.minPval # limits
        if elang<self.minPval: elang = self.minPval # limits
        if shang>self.maxPval: shang = self.maxPval # limits
        if elang>self.maxPval: elang = self.maxPval # limits
        elpos = [self.armLen[SH] * cos(shang), self.armLen[SH] * sin(shang)] # calculate shoulder x-y pos
        handpos = [elpos[X] + self.armLen[EL] * cos(shang+elang), elpos[Y] + self.armLen[EL] * sin(shang+elang)]
        shvel = self.angVel[SH] + (dataReceived[1]-dataReceived[0]) - (friction * self.angVel[SH])# update velocities based on incoming commands (accelerations) and friction
        elvel = self.angVel[EL] + (dataReceived[3]-dataReceived[2]) - (friction * self.angVel[EL])
        if self.anim:
            self.circle.center = self.targetPos
            self.line.set_data([0, elpos[0], handpos[0]], [0, elpos[1], handpos[1]]) # update line in figure
            self.ax.set_title('Time = %.1f ms, shoulder: pos=%.2f rad, vel=%.2f, acc=%.2f ; elbow: pos = %.2f rad, vel = %.2f, acc=%.2f' % (float(self.duration), shang, shvel, dataReceived[0] - (friction * shvel), elang, elvel, dataReceived[1] - (friction * elvel) ), fontsize=10)
            show(block=False)
            pause(0.0001) # pause so that the figure refreshes at every time step
        return [shang, elang, shvel, elvel, handpos[0], handpos[1]]

    #%% Reset arm variables so doesn't move between trials
    def resetArm(self, s, t):
        self.trial = self.trial + 1
        s.timeoflastreset = t
        if self.type == 'dummyArm':
            self.ang = list(self.startAng) # keeps track of shoulder and elbow angles
            self.angVel = [0,0] # keeps track of joint angular velocities
            self.motorCmd = [0,0,0,0] # motor commands to muscles
            self.error = 0 # error signal (eg. difference between )
            self.critic = 0 # critic signal (1=reward; -1=punishment)
            self.initArmMovement = self.initArmMovement + s.testTime
        

    def setPMdInput(self, s):
        for c in range(0, s.popnumbers[s.PMd], 2):
            try: 
                gid = s.popGidStart[s.PMd] + c # find gid of PMd 
                lid = s.gidDic[gid] # find local index corresponding to gid
                if (gid in s.targetPMdInputs[s.targetid]):  # if gid in PMdinputs for this target
                    #print "yes:",gid
                    s.cells[lid].interval=1000/s.maxPMdRate # set low interval (in ms as a function of rate)
                else: # if angle not in range -> low firing rate
                    s.cells[lid].interval=1000/s.minPMdRate # set high interval (in ms as a function of rate)
                    #print "no:",gid
            except:
                pass # local index corresponding to gid not found in this node

    #%% plot motor commands
    def RLcritic(self, t):
        if t > self.initArmMovement: # do not calculate critic signal in between trials
            # Calculate critic signal and activate RL (check synapses between nsloc->cells)
            RLsteps = int(self.RLinterval/self.interval)

            if len(self.errorAll) >= RLsteps:
                diff = self.error - mean(self.errorAll[-RLsteps:-1]) # difference between error at t and at t-(RLdt/dt) eg. t-50/5 = t-10steps
            else:
                diff = 0
            if diff < -self.minRLerror: # if error negative: LTP 
                self.critic = 1 # return critic signal to model.py so can update synaptic weights
            elif diff > self.minRLerror: # if error positive: LTD
                self.critic = -1
            else: # if difference not significant: no weight change
                self.critic = 0
        else: # if 
            self.critic = 0
        return self.critic

    #%% plot joint angles
    def plotTraj(self, filename):
        fig = figure() 
        l = 1.1*sum(self.armLen)
        ax = fig.add_subplot(111, autoscale_on=False, xlim=(-l/2, +l), ylim=(-l/2, +l)) # create subplot
        posX, posY = zip(*[self.angles2pos([x[SH],x[EL]], self.armLen) for x in self.angAll])
        ax.plot(posX, posY, 'r')
        targ = Circle((self.targetPos),0.04, color='g', fill=False) # target
        ax.add_artist(targ)
        ax.grid()
        ax.set_title('X-Y Hand trajectory')
        xlabel('x')
        ylabel('y')
        print 'saving '+filename
        fig.savefig(filename)

    #%% plot joint angles
    def plotAngs(self):
        fig = figure() 
        ax = fig.add_subplot(111) # create subplot
        sh = [x[SH] for x in self.angAll]
        el = [x[EL] for x in self.angAll]
        ax.plot(sh, 'r', label='shoulder')
        ax.plot(el, 'b', label='elbow')
        shTarg = self.pos2angles(self.targetPos, self.armLen)[0]
        elTarg = self.pos2angles(self.targetPos, self.armLen)[1]
        ax.plot(range(0,len(sh)), [shTarg] * len(sh), 'r:', label='sh target')
        ax.plot(range(0,len(el)), [elTarg] * len(el), 'b:', label='el target')
        ax.set_title('Joint angles')
        xlabel('time')
        ylabel('angle')
        ax.legend()

    #%% plot motor commands
    def plotMotorCmds(self):
        fig = figure() 
        ax = fig.add_subplot(111) # create subplot
        shext = [x[SH_EXT] for x in self.motorCmdAll]
        elext = [x[EL_EXT] for x in self.motorCmdAll]
        shflex = [x[SH_FLEX] for x in self.motorCmdAll]
        elflex = [x[EL_FLEX] for x in self.motorCmdAll]
        ax.plot(shext, 'r', label='sh ext')
        ax.plot(shflex, 'r:', label='sh flex')
        ax.plot(elext, 'b', label='el ext')
        ax.plot(elflex, 'b:', label='el flex')
        ax.set_title('Motor commands')
        xlabel('time')
        ylabel('motor command')
        ax.legend()

    #%% plot RL critic signal and error
    def plotRL(self):
        fig = figure() 
        ax = fig.add_subplot(111) # create subplot
        ax.plot(self.errorAll, 'r', label='error')
        ax.plot((array(self.criticAll)+1.0) * max(self.errorAll) / 2.0, 'b', label='RL critic')
        ax.set_title('RL critic and error')
        xlabel('time')
        ylabel('Error (m) / RL signal')
        ax.legend()

    ################################
    ### SETUP
    ################################
    def setup(self, s):#, nduration, loopstep, RLinterval, pc, scale, popnumbers, p): 
        self.duration = s.duration#/1000.0 # duration in msec
        self.interval = s.loopstep#/1000.0 # interval between arm updates in ,sec       
        self.RLinterval = s.RLinterval # interval between RL updates in msec
        self.minRLerror = s.minRLerror # minimum error change for RL (m)
        self.armLen = s.armLen # elbow - shoulder from MSM;radioulnar - elbow from MSM;  
        self.handPos = [0,0] # keeps track of hand (end-effector) x,y position
        self.handPosAll = [] # list with all handPos
        self.handVel = [0,0] # keeps track of hand (end-effector) x,y velocity
        self.handVelAll = [] # list with all handVel
        self.startAng = s.startAng # starting shoulder and elbow angles (rad) = natural rest position
        self.ang = list(self.startAng) # keeps track of shoulder and elbow angles
        self.angAll = [] # list with all ang
        self.angVel = [0,0] # keeps track of joint angular velocities
        self.angVelAll = [] # list with all angVel
        self.motorCmd = [0,0,0,0] # motor commands to muscles
        self.motorCmdAll = [] # list with all motorCmd
        self.targetDist = s.targetDist # target distance from center (15 cm)
        #self.targetid = 0 # target id (eg. 0=right, 1=left, 2=top, 3=bottom)
        self.targetidAll = [] # list with all targetid
        self.targetPos = self.setTargetByID(s.targetid, self.startAng, self.targetDist, self.armLen) # set the target location based on target id
        self.error = 0 # error signal (eg. difference between )
        self.errorAll = [] # list with all error
        self.critic = 0 # critic signal (1=reward; -1=punishment)
        self.criticAll = [] # list with all critic
        self.randDur = 0 # initialize explor movs duration
        self.initArmMovement = int(s.initArmMovement) # start arm movement after x msec
        self.trial = 0 # trial number

        # motor command encoding
        self.vec = h.Vector()
        self.cmdmaxrate = s.cmdmaxrate # maximum spikes for motor command (normalizing value)
        self.cmdtimewin = s.cmdtimewin # spike time window for shoulder motor command (ms)


        # proprioceptive encoding
        self.pStart = s.pStart # proprioceptive encoding first cell
        self.numPcells = s.numPcells # number of proprioceptive cells to encode shoulder and elbow angles
        self.prange = zeros((self.numPcells,2)) # range of values encoded by each P cell (proprioceptive tuning curve)
        self.minPval = s.minPval # min angle to encode
        self.maxPval = s.maxPval # max angle to encode
        self.minPrate = s.minPrate # firing rate when angle not within range
        self.maxPrate = s.maxPrate # firing rate when angle within range
        angInterval = (self.maxPval - self.minPval) / (self.numPcells - 2) * 2 # angle interval (times 2 because shoulder+elbow)
        currentPval = self.minPval
        for c in range(0,self.numPcells,2):
            self.prange[c,0] = currentPval # shoulder lower range
            self.prange[c,1] = currentPval + angInterval # shoulder higher range
            self.prange[c+1,0] = currentPval # elbow lower range
            self.prange[c+1,1] = currentPval + angInterval # elbow higher range
            currentPval += angInterval


        # initialize dummy or musculoskeletal arm 
        if s.rank == 0: 
            if self.type == 'dummyArm': 
                self.setupDummyArm() # setup dummyArm (eg. graph animation)
            elif self.type == 'musculoskeletal':  
                damping = 1 # damping of muscles (.osim parameter)
                shExtGain = 2  # gain factor to multiply force of shoulder extensor muscles (.osim parameter)
                shFlexGain = 1 # gain factor to multiply force of shoulder flexor muscles (.osim parameter)
                elExtGain = 1 # gain factor to multiply force of elbow extensor muscles (.osim parameter)
                elFlexGain = 0.8 # gain factor to multiply force of elbox flexor muscles (.osim parameter)
                # call function to initialize virtual arm params and run virtual arm C++ executable
                arminterface.setup(self.duration/1000.0, self.interval, self.startAng[SH], self.startAng[EL], self.targetPos[X], self.targetPos[Y], damping, shExtGain, shFlexGain, elExtGain, elFlexGain)
        
        # set PMd inputs
        if s.PMdinput == 'targetSplit': 
            self.setPMdInput(s) # set PMd inputs

    ################################          
    ### RUN     
    ################################
    def run(self, t, s): #pc, cells, gidVec, gidDic, cellsperhost=[], hostspikevecs=[]): 

        # Append to list the the value of relevant variables for this time step (only worker0)
        if s.rank == 0:     
            self.handPosAll.append(list(self.handPos)) # list with all handPos
            self.handVelAll.append(list(self.handVel))  # list with all handVel
            self.angAll.append(list(self.ang)) # list with all ang
            self.angVelAll.append(list(self.angVel)) # list with all angVel
            self.motorCmdAll.append(list(self.motorCmd)) # list with all motorCmd
            self.targetidAll.append(self.targetid) # list with all targetid
            self.errorAll.append(self.error) # error signal (eg. difference between )
            self.criticAll.append(self.critic) # critic signal (1=reward; -1=punishment)

        ############################
        # dummyArm or musculoskeletal: gather spikes for motor command
        ############################
        if self.type == 'dummyArm' or self.type == 'musculoskeletal': 
            ## Exploratory movements
            if s.explorMovs and t-s.timeoflastexplor >= self.randDur: # if time to update exploratory movement
                seed(s.id32('%d'%(int(t)+s.randseed))) # init seed
                self.randMus = int(uniform(0,s.nMuscles)) # select random muscle group
                self.randMul = uniform(s.explorMovsFactor/5,s.explorMovsFactor) # select random multiplier
                self.randDur = uniform(s.explorMovsDur/5, s.explorMovsDur) # select random duration
                s.timeoflastexplor = t

                if s.explorMovs == 1: # add random noise to EDSC+IDSC population
                    IDSCgids = array(s.motorCmdCellRange[self.randMus]) - int(s.popGidStart[s.EDSC]) + int(s.popGidStart[s.IDSC]) 
                    IDSCgidsInverse = [j for j in range(s.popGidStart[s.IDSC],s.popGidEnd[s.IDSC]+1) if j not in list(IDSCgids)]
                    for (i,x) in enumerate(s.backgroundsources): # for all background input netstims
                        if s.backgroundgid[i] in s.motorCmdCellRange[self.randMus] or s.backgroundgid[i] in list(IDSCgids): # if connected to chosen muscle cells
                            #print s.backgroundgid[i]
                            x.interval = (self.randMul*s.backgroundrateExplor)**-1*1e3  # increase firing
                        elif s.backgroundgid[i] in [y for z in range(s.nMuscles) if z != self.randMus for y in s.motorCmdCellRange[z]+IDSCgidsInverse]:
                            x.interval = s.backgroundrateMin**-1*1e3 # otherwise set to minimum
                    #if s.rank==0: print 'exploratory movement, randMus',self.randMus,' strength:',self.randMul,' duration:', self.randDur
                    #print 'IDSC cells:', IDSCgids   
                    #print 'IDSC inverse cells:', IDSCgidsInverse

                
                elif s.explorMovs == 2: # add random noise to EB5 population
                    self.targetCells = range(s.popGidStart[s.EB5], s.popGidEnd[s.EB5]+1) # EB5 cell gids
                    self.randNumCells = randint(1, int(s.explorCellsFraction*len(self.targetCells))) # num of cells to stimumales
                    self.randCells = sample(self.targetCells, int(self.randNumCells)) # select random gids
                    for (i,x) in enumerate(s.backgroundsources):  # for all input background netstims
                        if s.backgroundgid[i] in self.randCells:  # if connectd to selected random cells
                            x.interval = s.backgroundrateExplor**-1*1e3  # increase input
                        elif s.backgroundgid[i] in [y for y in self.targetCells if y not in self.randCells]:  # for the rest of cells
                            x.interval = s.backgroundrateMin**-1*1e3  # set to normal level
                    #if s.rank==0: 
                    #print 'Nodes:', s.rank,' - exploratory movement, numcells:',self.randNumCells,' strength:',self.randMul,' duration:', self.randDur, 'cells:', self.randCells    


            ## Reset arm and set target after every trial -start from center etc
            if s.trialReset and t-s.timeoflastreset > s.testTime: 
                self.resetArm(s, t)
                s.targetid = s.trialTargets[self.trial] # set target based on trial number
                self.targetPos = self.setTargetByID(s.targetid, self.startAng, self.targetDist, self.armLen) 
                if s.PMdinput == 'targetSplit': self.setPMdInput(s)


            ## Only move after initial period - avoids initial transitory spiking period (NSLOC sync spikes), and allows for variables with history to clear
            # can be justified as preparatory period (eg. watiing for go cue)
            if t > self.initArmMovement:
                ## Gather spikes #### from all vectors to then calculate motor command 
                for i in range(s.nMuscles):
                    cmdVecs = [array(s.hostspikevecs[x]) for x in range(s.cellsperhost) if (s.gidVec[x] in s.motorCmdCellRange[i])]
                    self.motorCmd[i] = sum([len(x[(x < t) * (x > t-self.cmdtimewin)]) for x in cmdVecs])
                    s.pc.allreduce(self.vec.from_python([self.motorCmd[i]]), 1) # sum
                    self.motorCmd[i] = self.vec.to_python()[0]       
            # else:
            #     for i in range(s.nMuscles): # stimulate all muscles equivalently so arm doesnt move
            #         self.motorCmd[i] = 0.2 * self.cmdmaxrate

                ## Calculate final motor command 
                if s.rank==0:  
                    self.motorCmd = [x / self.cmdmaxrate for x in self.motorCmd]  # normalize motor command 
                    if s.antagInh: # antagonist inhibition
                        if self.motorCmd[SH_EXT] > self.motorCmd[SH_FLEX]: # sh ext > sh flex
                            self.motorCmd[SH_FLEX] =  self.motorCmd[SH_FLEX]**2 / self.motorCmd[SH_EXT] / s.antagInh
                        elif self.motorCmd[SH_EXT] < self.motorCmd[SH_FLEX]: # sh flex > sh ext
                            self.motorCmd[SH_EXT] = self.motorCmd[SH_EXT]**2 / self.motorCmd[SH_FLEX] / s.antagInh
                        if self.motorCmd[EL_EXT] > self.motorCmd[EL_FLEX]: # el ext > el flex
                            self.motorCmd[EL_FLEX] = self.motorCmd[EL_FLEX]**2 / self.motorCmd[EL_EXT] / s.antagInh
                        elif self.motorCmd[EL_EXT] < self.motorCmd[EL_FLEX]: # el ext > el flex
                            self.motorCmd[EL_EXT] = self.motorCmd[EL_EXT]**2 / self.motorCmd[EL_FLEX] / s.antagInh
             

        ############################
        # ALL arms: Send motor command to virtual arm; receive new position; update proprioceptive population (ASC)
        ############################
        # Worker 0 sends motor commands and receives data from virtual arm
        #print "t=%f , self.initArmMovement=%f"%(t, self.initArmMovement)
        if s.rank == 0:
            if self.type == 'musculoskeletal': # MUSCULOSKELETAL
                try:
                    dataReceived = arminterface.sendAndReceiveDataPackets(t, self.interval, self.motorCmd[0], self.motorCmd[1], self.motorCmd[2], self.motorCmd[3])
                except:
                    dataReceived = [self.ang[SH], self.ang[EL]]
                if not dataReceived or dataReceived==[-3,3]:  # if error receiving packet
                    dataReceived = [self.ang[SH], self.ang[EL]]  # use previous packet
                    print 'Missed packet at t=%.2f',t
            elif self.type == 'dummyArm': # DUMMYARM
                dataReceived = self.runDummyArm(self.motorCmd) # run dummyArm
            elif self.type == 'randomOutput': # RANDOMOUTPUT
                dataReceived = [0,0] 
                dataReceived[0] = uniform(self.minPval, self.maxPval) # generate 2 random values  
                dataReceived[1] = uniform(self.minPval, self.maxPval)  
            # broadcast dataReceived  to other workers.   
            n = s.pc.broadcast(self.vec.from_python(dataReceived), 0) # convert python list to hoc vector for broadcast data received from arm
        else: # other workers
            n = s.pc.broadcast(self.vec, 0)  # receive shoulder and elbow angles from worker0 so can compare with cells in this worker
            dataReceived = self.vec.to_python()  
        if self.type == 'musculoskeletal':
            [self.ang[SH], self.ang[EL]] = dataReceived
            self.handPos = self.angles2pos(self.ang, self.armLen) 
            self.angVel[SH] = self.angVel[EL] = 0 
        else:
            [self.ang[SH], self.ang[EL], self.angVel[SH], self.angVel[EL], self.handPos[SH], self.handPos[EL]] = dataReceived # map data received to shoulder and elbow angles
        #[self.ang[SH], self.ang[EL], self.angVel[SH], self.angVel[EL]] = dataReceived # map data received to shoulder and elbow angles      
        
        #### Update proprio pop ASC
        for c in range(0, self.numPcells, 2):
            try: 
                id = s.gidDic[self.pStart + c] # find local index corresponding to gid
                if (self.ang[SH] >= self.prange[c,0] and self.ang[SH] < self.prange[c,1]):  # in angle in range -> high firing rate
                    s.cells[id].interval=1000/self.maxPrate # interval in ms as a function of rate
                else: # if angle not in range -> low firing rate
                    s.cells[id].interval=1000/self.minPrate # interval in ms as a function of rate
            except:
                pass # local index corresponding to gid not found in this node
            try: 
                id = s.gidDic[self.pStart + c + 1] # find local index corresponding to gid
                if (self.ang[EL] >= self.prange[c+1,0] and self.ang[EL] < self.prange[c+1,1]):  # in angle in range -> high firing rate
                    s.cells[id].interval=1000/self.maxPrate # interval in ms as a function of rate
                else: # if angle not in range -> low firing rate
                    s.cells[id].interval=1000/self.minPrate # interval in ms as a function of rate 
            except:
                pass # local index corresponding to gid not found in this node


        #### Calculate error between hand and target for interval between RL updates 
        if s.rank == 0 and self.initArmMovement: # do not update between trials
            #print 't=%.2f, xpos=%.2f'%(t,self.targetPos[X])
            self.error = sqrt((self.handPos[X] - self.targetPos[X])**2 + (self.handPos[Y] - self.targetPos[Y])**2)
            
        return self.critic


    ################################
    ### CLOSE
    ################################
    def close(self, s):             
        if self.type == 'randomOutput':
            print('\nClosing random output virtual arm...')

        if self.type == 'dummyArm':
            if s.explorMovs == 1: # remove explor movs related noise to cells
                for imus in range(s.nMuscles):
                    IDSCgids = array(s.motorCmdCellRange[self.randMus]) - int(s.popGidStart[s.EDSC]) + int(s.popGidStart[s.IDSC]) 
                    for (i,x) in enumerate(s.backgroundsources):
                        if s.backgroundgid[i] in s.motorCmdCellRange[imus]+list(IDSCgids):
                            x.interval = 0.0001**-1*1e3
                            x.noise = s.backgroundnoise # Fractional noise in timing
            elif s.explorMovs == 2: # remove explor movs related noise to cells
                for (i,x) in enumerate(s.backgroundsources):  # for all input background netstims
                    if s.backgroundgid[i] in self.targetCells:  # for the rest of cells
                        x.interval = s.backgroundrate**-1*1e3  # set to normal level

            if s.trialReset:
                s.timeoflastreset = 0
                self.initArmMovement = int(s.initArmMovement)

            if s.rank == 0:
                print('\nClosing dummy virtual arm ...') 

                if self.anim:
                    ioff() # turn interactive mode off
                    close(self.fig) # close arm animation graph 
                if self.graphs: # plot graphs
                    self.plotTraj(s.filename)
                    #self.plotAngs()
                    #self.plotMotorCmds()
                    #self.plotRL()
        
        if self.type == 'musculoskeletal':
            if s.explorMovs == 1: # remove explor movs related noise to cells
                for imus in range(s.nMuscles):
                    IDSCgids = array(s.motorCmdCellRange[self.randMus]) - int(s.popGidStart[s.EDSC]) + int(s.popGidStart[s.IDSC]) 
                    for (i,x) in enumerate(s.backgroundsources):
                        if s.backgroundgid[i] in s.motorCmdCellRange[imus]+list(IDSCgids):
                            x.interval = 0.0001**-1*1e3
                            x.noise = s.backgroundnoise # Fractional noise in timing
            elif s.explorMovs == 2: # remove explor movs related noise to cells
                for (i,x) in enumerate(s.backgroundsources):  # for all input background netstims
                    if s.backgroundgid[i] in self.targetCells:  # for the rest of cells
                        x.interval = s.backgroundrate**-1*1e3  # set to normal level

            if s.trialReset:
                s.timeoflastreset = 0
                self.initArmMovement = int(s.initArmMovement)

            if s.rank == 0:
                print('\nClosing dummy virtual arm ...') 
                arminterface.closeSavePlot(self.duration/1000.0, self.interval)

                if self.graphs: # plot graphs
                    self.plotTraj(s.filename)
                    #self.plotAngs()
                    self.plotMotorCmds()
                    self.plotRL()