GSoC Open Source Brain: Cortical Connections

Cortical Connections

In the same vein that the post before this one we will show here how to construct the connections between the cortical layers. In order to do so we will construct a function that works in general for any arbitrary connectivity, we describe in the following its structure. First, as in the thalamo-cortical connectivity, we have again the same structure of a function that loops over the target population extracting the relevant parameters that characterize these neurons. Furthermore we have another function that loops over the source population creating the corresponding tuples for the connection list. Is in this last function where the particular connectivity rule is implemented.

In the particular case of the Troyer model the connectivity between the cortical cells is determined by the correlation between the receptive fields of the neurons, the receptive fields here being Gabor functions. In more detail the neurons whose receptive fields are more correlated will be the ones more likely to have excitatory connections between them. On the other hand the ones whose receptive fields are less correlated will be more likely to receive inhibitory connections. In this post we show two schemes that accomplish this connectivity. The first one uses the fact the parameters of the receptive field to calculate a connectivity and the second one uses the receptive fields directly to calculate the correlations. We present the determining functions in the stated order down here.

Now we present the function that creates the connectivity for a given neuron par. The circular distance between the orientation and phases are calculated as a proxy to estimate how similar the receptive fields of the neurons are. After that, the distance between them is weighted and normalized with a normal function in order to obtain a value that we can interpret as a probability value. Finally in order to calculate the connectivity we sample n_pick times with the given probability value to see hwo strong a particular connection should be.

 
def cortical_to_cortical_connection(target_neuron_index, connections, source_population, n_pick, g, delay, source_orientations,
                                    source_phases, orientation_sigma, phase_sigma, target_neuron_orientation,
                                    target_neuron_phase, target_type):
    """
    Creates the connections from the source population to the target neuron

    """
    for source_neuron in source_population:
        # Extract index, orientation and phase of the target
        source_neuron_index = source_population.id_to_index(source_neuron)
        source_neuron_orientation = source_orientations[source_neuron_index]
        source_neuron_phase = source_phases[source_neuron_index]

        # Now calculate phase and orientation distances
        or_distance = circular_dist(target_neuron_orientation, source_neuron_orientation, 180)

        if target_type:
            phase_distance = circular_dist(target_neuron_phase, source_neuron_phase, 360)
        else:
            phase_distance = 180 - circular_dist(target_neuron_phase, source_neuron_phase, 360)

        # Now calculate the gaussian function
        or_gauss = normal_function(or_distance, mean=0, sigma=orientation_sigma)
        phase_gauss = normal_function(phase_distance, mean=0, sigma=phase_sigma)

        # Now normalize by guassian in zero
        or_gauss = or_gauss / normal_function(0, mean=0, sigma=orientation_sigma)
        phase_gauss = phase_gauss / normal_function(0, mean=0, sigma=phase_sigma)

        # Probability is the product
        probability = or_gauss * phase_gauss
        probability = np.sum(np.random.rand(n_pick) < probability)  # Samples
        synaptic_weight = (g / n_pick) * probability


        if synaptic_weight > 0:
                    connections.append((source_neuron_index, target_neuron_index, synaptic_weight, delay))

    return connections
      
    

Note that the overall strength is weighted by the conductivity value g that is passed as an argument. Furthermore a delay that is also passed as an argument is added to the list as the last element of the tuple.

Secondly we present the full correlation scheme. In this scheme we utilize the kernels directly to calculate the spatial correlation between them. In particular after we have flattened our kernels Z to have a series instead of a matrix we use the function perasonr from scipy.stats to calculate the correlation. Again as in the case above we use this probability to sample n_pick times and then calculate the relative connectivity strength with this.

 
def cortical_to_cortical_connection_corr(target_neuron_index, connections, source_population, n_pick, g, delay,
                                    source_orientations, source_phases, target_neuron_orientation, target_neuron_phase,
                                    Z1, lx, dx, ly, dy, sigma, gamma, w, target_type):
    """
    Creates the connections from the source population to the target neuron

    """
    for source_neuron in source_population:
        # Extract index, orientation and phase of the target
        x_source, y_source = source_neuron.position[0:2]
        source_neuron_index = source_population.id_to_index(source_neuron)
        source_neuron_orientation = source_orientations[source_neuron_index]
        source_neuron_phase = source_phases[source_neuron_index]

        Z2 = gabor_kernel(lx, dx, ly, dy, sigma, gamma, source_neuron_phase, w, source_neuron_orientation,
                          x_source, y_source)


        if target_type:
            probability = pearsonr(Z1.flat, Z2.flat)[0]
        else:
            probability = (-1) * pearsonr(Z1.flat, Z2.flat)[0]

        probability = np.sum(np.random.rand(n_pick) < probability)  # Samples
        synaptic_weight = (g / n_pick) * probability

        if synaptic_weight > 0:
                    connections.append((source_neuron_index, target_neuron_index, synaptic_weight, delay))

    return connections
    

Note that the overall strength is weighted by the conductivity value g that is passed as an argument. Furthermore a delay that is also passed as an argument is added to the list as the last element of the tuple.

We now show how a plot that illustrates how the probabilities change when the parameters that determined the gabor function are changed for each scheme.

In the figure above we have int he upper part how the probability for the first scheme a neuron with phase 0 and orientation 0 change as we vary the phase (left) and orientation (right). In the two graphs bellow we have the same for the second scheme we presented

GSoC Open Source Brain: Thalamo-Cortical Connections

Thalamo-cortical connections

In this post I will show how to build arbitrary custom connections in PyNN. We will illustrate the general technique in the particular case of the Troyer model. In the Troyer model the connections from the LGN to the cortex are determined with a gabor-profile therefore I am going to describe the required functions to achieve such an aim.

In the PyNN documentation we find that one of the ways of implementing arbitrary connectivity patterns is to use the FromListConnector utility. In this format we have to construct a list of tuples with a tuple for each connection. In each tuple we need to include the index of the source neuron (the neuron from which the synapse originates), the index of the target neuron (the neuron into which the synapse terminates), the weight and the delay. For example (0, 1, 5, 0.1) would indicate that we have a connection from the neuron 0 to the neuron 1 with a synaptic weight of 5 and a delay of 0.1.

In the light of the explanation above we need to construct a function that is able to construct a list with the appropriate weights given a target and a source populations. In order to start moving towards this goal we will first write a function that connects a given neuron in the target population to all the neurons in the source population. We first present the function here bellow and we will explain it later:

 
def lgn_to_cortical_connection(cortical_neuron_index, connections, lgn_neurons, n_pick, g, delay, polarity, sigma,
                               gamma, phi, w, theta, x_cortical, y_cortical):
    """
    Creates connections from the LGN to the cortex with a Gabor profile.

    This function adds all the connections from the LGN to the cortical cell with index = cortical_neuron_index. It
    requires as parameters the cortical_neruon_index, the current list of connections, the lgn population and also
    the parameters of the Gabor function.

    Parameters
    ----
    cortical_neuron_index : the neuron in the cortex -target- that we are going to connect to
    connections: the list with the connections to which we will append the new connnections
    lgn_neurons: the source population
    n_pick: How many times we will sample per neuron
    g: how strong is the connection per neuron
    delay: the time it takes for the action potential to arrive to the target neuron from the source neuron
    polarity: Whether we are connection from on cells or off cells
    sigma: Controls the decay of the exponential term
    gamma: x:y proportionality factor, elongates the pattern
    phi: Phase of the overall pattern
    w: Frequency of the pattern
    theta: Rotates the whole pattern by the angle theta
    x_cortical, y_cortical : The spatial coordinate of the cortical neuron

    """

    for lgn_neuron in lgn_neurons:
            # Extract position
            x, y = lgn_neuron.position[0:2]
            # Calculate the gabbor probability
            probability = polarity * gabor_probability(x, y, sigma, gamma, phi, w, theta, x_cortical, y_cortical)
            probability = np.sum(np.random.rand(n_pick) < probability)  # Samples

            synaptic_weight = (g / n_pick) * probability
            lgn_neuron_index = lgn_neurons.id_to_index(lgn_neuron)

            # The format of the connector list should be pre_neuron, post_neuron, w, tau_delay
            if synaptic_weight > 0:
                connections.append((lgn_neuron_index, cortical_neuron_index, synaptic_weight, delay))

    

The first thing to note from the function above are its arguments. It contains the source population and the particular target neuron that we want to connect to. It also contains all the connectivity and gabor-function related parameters. In the body of the function we have one loop over the whole source population that decides whether we add a connection from a particular cell or not. In order to decide if we add a connection we have to determine the probability from the gabor function. Once we have this we sample n_pick times and add a weighted synaptic weight accordingly to the list for each neuron.

In the function above we have the values of the gabor function passed as arguments. However, the values of each gabor function depend on the nature of the cell of the target population. In the light of this we will construct another function that loops over the target population and extracts the appropriate gabor values for each function in this population. We again present the function and then explain it:

def create_lgn_to_cortical(lgn_population, cortical_population, polarity,  n_pick, g, delay,  sigma, gamma, phases,
                           w, orientations):
    """
    Creates the connection from the lgn population to the cortical population with a gabor profile. It also extracts
    the corresponding gabor parameters that are needed in order to determine the connectivity.
    """

    print 'Creating connection from ' + lgn_population.label + ' to ' + cortical_population.label

    # Initialize connections
    connections = []

    for cortical_neuron in cortical_population:
        # Set the parameters
        x_cortical, y_cortical = cortical_neuron.position[0:2]
        cortical_neuron_index = cortical_population.id_to_index(cortical_neuron)
        theta = orientations[cortical_neuron_index]
        phi = phases[cortical_neuron_index]

        # Create the connections from lgn to cortical_neuron
        #lgn_to_cortical_connection(cortical_neuron_index, connections, lgn_population, n_pick, g, polarity, sigma,
        #gamma, phi, w, theta, x_cortical, y_cortical)

        lgn_to_cortical_connection(cortical_neuron_index, connections, lgn_population, n_pick, g, delay, polarity, sigma,
                                   gamma, phi, w, theta, 0, 0)

    return connections
    

This function requires as arguments the source and target populations as well as the necessary parameters that characterize each cell connectivity: orientation and phase. In the body of the function we have a loop over the cortical population that extracts the relevant parameters -position, orientation and phase- and then calls the function that we already describe previously in order to create the connectivity from the source population to the cell in place.

So now we have the necessary functions to construct a list. Now, we can use FromListConnector to transform the list into a connector. And the use this to define a Projection. We define both the excitatory and inhibitory connections. We abstract this complete set into the following function:

def create_thalamocortical_connection(source, target, polarity, n_pick, g, delay, sigma, gamma, w, phases, orientations, simulator):
    """
    Creates a connection from a layer in the thalamus to a layer in the cortex through the mechanism of Gabor sampling
    """

    # Produce a list with the connections
    connections_list = create_lgn_to_cortical(source, target, polarity, n_pick, g, delay, sigma, gamma, phases, w, orientations)

    # Transform it into a connector
    connector = simulator.FromListConnector(connections_list, column_names=["weight", "delay"])

    # Create the excitatory and inhibitory projections
    simulator.Projection(source, target, connector, receptor_type='excitatory')
    simulator.Projection(source, target, connector, receptor_type='inhibitory')      
    

With this we can create in general connections from one target population to the other. We can even change change the gabor function for whatever we want if we want to experiment with other connectivity patterns. Finally we present down here an example of a sampling from a Gabor function with the aglorithm we just constructed:

So in the image we show in the left the sampling from the ideal Gabor function in the right.

GSoC Open Source Brain: Arbitrary Spike-trains in PyNN

Arbitrary Spike-trains in PyNN

In this example we are going to create a population of cells with arbitrary spike trains. We will load the spike train from a file where they are stored as a list of arrays with the times at which they occurred. In order to so we are going to use the SpikeSourceArray class model of PyNN

First we start importing PyNN in with nest and all the other required libraries. Furthermore we start the simulators and given some general parameters. We assume that we already have produced spikes for the cells with 0.50 contrast:

      import numpy as np
      import matplotlib.pyplot as plt
      import cPickle
      import pyNN.nest as simulator

      contrast = 0.50
      Nside_lgn = 30
      Ncell_lgn = Nside_lgn * Nside_lgn
      N_lgn_layers = 4
      t = 1000  # ms

      simulator.setup(timestep=0.1, min_delay=0.1, max_delay=5.0)
    

So we are going to suppose that we have our data stored in './data'. The spike-trains are lists as long as the cell population that contain for each element an array with the times at which the spikes occurred for that particular neuron. In order to load them we will use the following code

      directory = './data/'
      format = '.cpickle'

      spikes_on = []
      spikes_off = []

      for layer in xrange(N_lgn_layers):

      #  Layer 1
      layer = '_layer' + str(layer)

      polarity = '_on'
      contrast_mark = str(contrast)
      mark = '_spike_train'
      spikes_filename = directory + contrast_mark + mark + polarity + layer + format
      f2 = open(spikes_filename, 'rb')
      spikes_on.append(cPickle.load(f2))
      f2.close()

      polarity = '_off'
      contrast_mark = str(contrast)
      mark = '_spike_train'
      spikes_filename = directory + contrast_mark + mark + polarity + layer + format
      f2 = open(spikes_filename, 'rb')
      spikes_off.append(cPickle.load(f2))
      f2.close()

    

Now this is the crucial part. If we want to utilize the SpikeSourceArray model for a cell in PyNN we can define a function that pass the spike-train for each cell in the population. In order to so we use the following code:

      def spike_times(simulator, layer, spikes_file):
         return [simulator.Sequence(x) for x in spikes_file[layer]]
    

Note that we have to change every spike-train array to a sequence before using it as a spike-train. After defining this function we can create the LGN models:

      # Cells models for the LGN spikes (SpikeSourceArray)
      lgn_spikes_on_models = []
      lgn_spikes_off_models = []


      for layer in xrange(N_lgn_layers):
         model = simulator.SpikeSourceArray(spike_times=spike_times(simulator, layer, spikes_on))
         lgn_spikes_on_models.append(model)
         model = simulator.SpikeSourceArray(spike_times=spike_times(simulator, layer, spikes_off))
         lgn_spikes_off_models.append(model)
    

Now that we have the corresponding model for the cells we can create the populations in the usual way:

 # LGN Popluations

 lgn_on_populations = []
 lgn_off_populations = []

 for layer in xrange(N_lgn_layers):
 population = simulator.Population(Ncell_lgn, lgn_spikes_on_models[layer], label='LGN_on_layer_' + str(layer))
 lgn_on_populations.append(population)
 population = simulator.Population(Ncell_lgn, lgn_spikes_off_models[layer], label='LGN_off_layer_' + str(layer))
 lgn_off_populations.append(population)
      

In order to analyze the spike-trains patterns for each population we need to declare a recorder for each population:

 
      layer = 0  # We declare here the layer of our interest 

      population_on = lgn_on_populations[layer]
      population_off = lgn_off_populations[layer]

      population_on.record('spikes')
      population_off.record('spikes')
    

Note here that we can chose the layer of our interest by modifying the value of the layer variable. Finally we run the model with the usual instructions and extract the spikes:

      #############################
      # Run model
      #############################

      simulator.run(t)  # Run the simulations for t ms
      simulator.end()

      #############################
      # Extract the data
      #############################
      data_on = population_on.get_data()  # Creates a Neo Block
      data_off = population_off.get_data()

      segment_on = data_on.segments[0]  # Takes the first segment
      segment_off = data_off.segments[0]
    

In order to visualize the spikes we use the following function:

 
      # Plot spike trains
      def plot_spiketrains(segment):
          """
          Plots the spikes of all the cells in the given segments
          """
          for spiketrain in segment.spiketrains:
              y = np.ones_like(spiketrain) * spiketrain.annotations['source_id']
              plt.plot(spiketrain, y, '*b')
              plt.ylabel('Neuron number')
              plt.xlabel('Spikes')
    

Here the spiketrain variable contains the spike-train for each cell, that is, an array with the times at which the action potentials happened for each cell. In order to tell them apart we assigned them the value of the cell id. Finally we can plot the spikes of the on and off cells with the following code:

      plt.subplot(2, 1, 1)
      plt.title('On cells ')
      plot_spiketrains(segment_on)

      plt.subplot(2, 1, 2)
      plt.title('Off cells ')
      plot_spiketrains(segment_off)

      plt.show()
    

We now show the plot produced by the code above. Note that the on and off cells are off-phase by 180.

GSoC Open Source Brain: Firing Rate Induced by a Sinus Grating

Firing Rate induced by a Sinus Grating

Now that we know how to do convolutions with our center-surround kernel we can chose any other kind of stimulus to carry this out. In the neuoscience of vision it is very common to use a sinus grating in a wide array of experimental setings so we are going to use it now. In short, in this post we are going to see the see what signal does a center-surround kernel produces when is convolved with a sinus grating.

Center-Surround Kernel

In order to do the convolution we are going to define the kernel in the usual way using a function that we have utilized from our work before:

      # First we define the size and resolution of the space in which the convolution is going to happen
      dx = 0.05
      dy = 0.05
      lx = 6.0  # In degrees
      ly = 6.0  # In degrees

      # Now we define the temporal parameters of the kernel
      dt_kernel = 5.0  # ms
      kernel_duration = 150  # ms
      kernel_size = int(kernel_duration / dt_kernel)

      #  Now the center surround parameters
      factor = 1  # Controls the overall size of the center-surround pattern
      sigma_center = 0.25 * factor  # Corresponds to 15'
      sigma_surround = 1 * factor  # Corresponds to 1 degree

      # Finally we create the kernel
      kernel_on = create_kernel(dx, lx, dy, ly, sigma_surround, sigma_center, dt_kernel, kernel_size)
    

Sinus Grating

Now we are going to construct our sinus grating. But first, we need to think on how long our stimulus is going to last which is a function of how long the we want to simulate the convolution and of the resolutions of the stimulus and the simulation:

      ## Now we define the temporal l parameters of the sinus grating
      dt_stimuli = 5.0  # ms

      # We also need to add how long do we want to convolve
      dt = 1.0  # Simulation resolution
      T_simulation = 1 * 10 ** 3.0 # ms
      T_simulation += int(kernel_size * dt_kernel)  # Add the size of the kernel
      Nt_simulation = int(T_simulation / dt)  # Number of simulation points
      N_stimuli = int(T_simulation / dt_stimuli)  # Number of stimuli points     
    

Finally we now present the parameters that determine the sinus grating. First the spatial frequency (K), followed by the spatial phase (Phi) and orientation (Theta). Furthermore we have also a parameter for the amplitude and the temporal frequency:

      # And now the spatial parameters of the sinus grating
      K = 0.8  # Cycles per degree
      Phi = 0  # Spatial phase 
      Theta = 0 # Orientation 
      A = 1 # Amplitude 
      # Temporal frequency of sine grating
      w = 3  # Hz
    

Now with all the spatial parameters in our possession we can call the function that produces the sine grating, we define it as the following function that we present below:

 
      stimuli = sine_grating(dx, lx, dy, ly, A, K, Phi, Theta, dt_stimuli, N_stimuli, w)
      
      def sine_grating(dx, Lx, dy, Ly, A, K, Phi, Theta, dt_stimuli, N_stimuli, w):
       '''
       Returns a sine grating stimuli
       '''
       Nx = int(Lx / dx)
       Ny = int(Ly / dy)

       # Transform to appropriate units
       K = K * 2 * np.pi # Transforms K to cycles per degree
       w = w / 1000.0 # Transforms w to kHz

       x = np.arange(-Lx/2, Lx/2, dx)
       y = np.arange(-Ly/2, Ly/2, dy)
       X, Y = np.meshgrid(x, y)
       Z = A * np.cos(K * X *cos(Theta) + K * Y * sin(Theta) - Phi)
       t = np.arange(0, N_stimuli * dt_stimuli, dt_stimuli)
       f_t = np.cos(w * 2 * np.pi *  t )

       stimuli = np.zeros((N_stimuli, Nx, Ny))

       for k, time_component in enumerate(f_t):
           stimuli[k, ...] = Z * time_component

       return stimuli

    

Convolution

Now that we have the stimulus and the kernel we can do the convolution, in order to do that we use again our functions and indexes that we use in the last post:

 
      ## Now we can do the convolution

      # First we define the necessary indexes to the convolution
      signal_indexes, delay_indexes, stimuli_indexes = create_standar_indexes(dt, dt_kernel, dt_stimuli, kernel_size, Nt_simulation)
      working_indexes, kernel_times = create_extra_indexes(kernel_size, Nt_simulation)

      # Now we calculate the signal
      signal = np.zeros(Nt_simulation)

      for index in signal_indexes:
          signal[index] = convolution(index, kernel_times, delay_indexes, stimuli_indexes, kernel_on, stimuli)
    

We can visualize signal with the following code:

 
      #Plot the signal 
      t = np.arange(kernel_size*dt_kernel, T_simulation, dt)
      plt.plot(t, signal[signal_indexes])
      plt.show()
    

We can see that the signal is also a sinus with a frequency that is consistent with the one from the sinus grating.